blob: 33182c0562b8302aa2464ebf6b80b1d6a3a02aa9 [file] [log] [blame]
 /*************************************************************** ****************************************************************** **** Gaussian Elimination with partial pivoting. **** **** This file contains the factorization routine SGEFA **** ****************************************************************** ****************************************************************/ #include "ge.h" static char SGEFASid[] = "@(#)sgefa.c 1.1 2/4/86"; int sgefa( a, ipvt ) struct FULL *a; int *ipvt; /* PURPOSE SGEFA factors a real matrix by gaussian elimination. REMARKS SGEFA is usually called by SGECO, but it can be called directly with a saving in time if rcond is not needed. (time for SGECO) = (1 + 9/n)*(time for SGEFA) . INPUT a A pointer to the FULL matrix structure. See the definition in ge.h. OUTPUT a A pointer to the FULL matrix structure containing an upper triangular matrix and the multipliers which were used to obtain it. The factorization can be written a = l*u where l is a product of permutation and unit lower triangular matrices and u is upper triangular. ipvt An integer vector (of length a->cd) of pivot indices. RETURNS = -1 Matrix is not square. = 0 Normal return value. = k if u(k,k) .eq. 0.0 . This is not an error condition for this subroutine, but it does indicate that sgesl or sgedi will divide by zero if called. Use rcond in sgeco for a reliable indication of singularity. ROUTINES blas ISAMAX */ { register int i, j; int isamax(), k, l, nm1, info, n; float t, *akk, *alk, *aij, *mik; /* Gaussian elimination with partial pivoting. */ if( a->cd != a->rd ) return( -1 ); n = a->cd; nm1 = n - 1; akk = a->pd; info = 0; /* Assume nothing will go wrong! */ if( n < 2 ) goto CLEAN_UP; /* Loop over Diagional */ for( k=0; kpd[k] + k; l = isamax( n-k, akk, 1 ) + k; *ipvt = l; /* Zero pivot implies this column already triangularized. */ alk = a->pd[k] + l; if( *alk == 0.0e0) { info = k; continue; } /* Interchange a(k,k) and a(l,k) if necessary. */ if( l != k ) { t = *alk; *alk = *akk; *akk = t; } /* Compute multipliers for this column. */ t = -1.0e0 / (*akk); for( i=k+1, mik = akk+1; ipd[j]+k+1, mik=akk+1; i