*DECK DGEFA
SUBROUTINE DGEFA (A, LDA, N, IPVT, INFO)
C***BEGIN PROLOGUE DGEFA
C***PURPOSE Factor a matrix using Gaussian elimination.
C***LIBRARY SLATEC (LINPACK)
C***CATEGORY D2A1
C***TYPE DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C)
C***KEYWORDS GENERAL MATRIX, LINEAR ALGEBRA, LINPACK,
C MATRIX FACTORIZATION
C***AUTHOR Moler, C. B., (U. of New Mexico)
C***DESCRIPTION
C
C DGEFA factors a double precision matrix by Gaussian elimination.
C
C DGEFA is usually called by DGECO, but it can be called
C directly with a saving in time if RCOND is not needed.
C (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .
C
C On Entry
C
C A DOUBLE PRECISION(LDA, N)
C the matrix to be factored.
C
C LDA INTEGER
C the leading dimension of the array A .
C
C N INTEGER
C the order of the matrix A .
C
C On Return
C
C A an upper triangular matrix and the multipliers
C which were used to obtain it.
C The factorization can be written A = L*U where
C L is a product of permutation and unit lower
C triangular matrices and U is upper triangular.
C
C IPVT INTEGER(N)
C an integer vector of pivot indices.
C
C INFO INTEGER
C = 0 normal value.
C = K if U(K,K) .EQ. 0.0 . This is not an error
C condition for this subroutine, but it does
C indicate that DGESL or DGEDI will divide by zero
C if called. Use RCOND in DGECO for a reliable
C indication of singularity.
C
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX
C***REVISION HISTORY (YYMMDD)
C 780814 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DGEFA
INTEGER LDA,N,IPVT(*),INFO
DOUBLE PRECISION A(LDA,*)
C
DOUBLE PRECISION T
INTEGER IDAMAX,J,K,KP1,L,NM1
C
C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
C***FIRST EXECUTABLE STATEMENT DGEFA
INFO = 0
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 70
DO 60 K = 1, NM1
KP1 = K + 1
C
C FIND L = PIVOT INDEX
C
L = IDAMAX(N-K+1,A(K,K),1) + K - 1
IPVT(K) = L
C
C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
IF (A(L,K) .EQ. 0.0D0) GO TO 40
C
C INTERCHANGE IF NECESSARY
C
IF (L .EQ. K) GO TO 10
T = A(L,K)
A(L,K) = A(K,K)
A(K,K) = T
10 CONTINUE
C
C COMPUTE MULTIPLIERS
C
T = -1.0D0/A(K,K)
CALL DSCAL(N-K,T,A(K+1,K),1)
C
C ROW ELIMINATION WITH COLUMN INDEXING
C
DO 30 J = KP1, N
T = A(L,J)
IF (L .EQ. K) GO TO 20
A(L,J) = A(K,J)
A(K,J) = T
20 CONTINUE
CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
30 CONTINUE
GO TO 50
40 CONTINUE
INFO = K
50 CONTINUE
60 CONTINUE
70 CONTINUE
IPVT(N) = N
IF (A(N,N) .EQ. 0.0D0) INFO = N
RETURN
END