blob: bd38c5a318138e757e05ba1af6e904854964925e [file] [log] [blame]
 /*************************************************************** ****************************************************************** **** Gaussian Elimination with partial pivoting. **** **** This file contains the solution routine SGESL **** ****************************************************************** ****************************************************************/ #include "ge.h" static char SGESLSid[] = "@(#)sgesl.c 1.1 2/4/86"; int sgesl( a, ipvt, b, job ) struct FULL *a; int *ipvt, job; float b[]; /* PURPOSE SGESL solves the real system a * x = b or trans(a) * x = b using the factors computed by SGECO or SGEFA. INPUT a A pointer to the FULL matrix structure containing the factored matrix. See the definition of FULL in ge.h. ipvt The pivot vector (of length a->cd) from SGECO or SGEFA. b The right hand side vector (of length a->cd). job = 0 to solve a*x = b , = nonzero to solve trans(a)*x = b where trans(a) is the transpose. OUTPUT b The solution vector x. REMARKS Error condition: A division by zero will occur if the input factor contains a zero on the diagonal. Technically this indicates singularity but it is often caused by improper arguments or improper setting of lda . It will not occur if the subroutines are called correctly and if sgeco has set rcond .gt. 0.0 or sgefa has set info .eq. 0 . */ { float t; float *akk, *mik, *uik, *bi; register int i, k; int l, n, nm1; n = a->cd; nm1 = n - 1; /* job = 0 , solve a * x = b. */ if( job == 0 ) { /* Forward elimination. Solve l*y = b. */ for( k=0; kpd[k] + k; /* akk points to a(k,k). */ /* Interchange b[k] and b[l] if necessary. */ l = *ipvt; t = b[l]; if( l != k ) { b[l] = b[k]; b[k] = t; } for( i=k+1, mik=akk+1; i=0; k-- ) { akk = a->pd[k] +k; b[k] /= (*akk); for( i=0, uik=a->pd[k]; ipd[k] + k; for( i=0, t=0.0, uik=a->pd[k], bi=b; i=0; k--, ipvt-- ) { for( i=k+1, t=0.0, mik=a->pd[k]+k+1, bi=b+k+1; i