share/odepack/fortran/dgefa.f

   1 *DECK DGEFA
   2       SUBROUTINE DGEFA (A, LDA, N, IPVT, INFO)
   3 C***BEGIN PROLOGUE  DGEFA
   4 C***PURPOSE  Factor a matrix using Gaussian elimination.
   5 C***CATEGORY  D2A1
   6 C***TYPE      DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C)
   7 C***KEYWORDS  GENERAL MATRIX, LINEAR ALGEBRA, LINPACK,
   8 C             MATRIX FACTORIZATION
   9 C***AUTHOR  Moler, C. B., (U. of New Mexico)
  10 C***DESCRIPTION
  11 C
  12 C     DGEFA factors a double precision matrix by Gaussian elimination.
  13 C
  14 C     DGEFA is usually called by DGECO, but it can be called
  15 C     directly with a saving in time if  RCOND  is not needed.
  16 C     (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .
  17 C
  18 C     On Entry
  19 C
  20 C        A       DOUBLE PRECISION(LDA, N)
  21 C                the matrix to be factored.
  22 C
  23 C        LDA     INTEGER
  24 C                the leading dimension of the array  A .
  25 C
  26 C        N       INTEGER
  27 C                the order of the matrix  A .
  28 C
  29 C     On Return
  30 C
  31 C        A       an upper triangular matrix and the multipliers
  32 C                which were used to obtain it.
  33 C                The factorization can be written  A = L*U  where
  34 C                L  is a product of permutation and unit lower
  35 C                triangular matrices and  U  is upper triangular.
  36 C
  37 C        IPVT    INTEGER(N)
  38 C                an integer vector of pivot indices.
  39 C
  40 C        INFO    INTEGER
  41 C                = 0  normal value.
  42 C                = K  if  U(K,K) .EQ. 0.0 .  This is not an error
  43 C                     condition for this subroutine, but it does
  44 C                     indicate that DGESL or DGEDI will divide by zero
  45 C                     if called.  Use  RCOND  in DGECO for a reliable
  46 C                     indication of singularity.
  47 C
  48 C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
  49 C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
  50 C***ROUTINES CALLED  DAXPY, DSCAL, IDAMAX
  51 C***REVISION HISTORY  (YYMMDD)
  52 C   780814  DATE WRITTEN
  53 C   890831  Modified array declarations.  (WRB)
  54 C   890831  REVISION DATE from Version 3.2
  55 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  56 C   900326  Removed duplicate information from DESCRIPTION section.
  57 C           (WRB)
  58 C   920501  Reformatted the REFERENCES section.  (WRB)
  59 C***END PROLOGUE  DGEFA
  60       INTEGER LDA,N,IPVT(*),INFO
  61       DOUBLE PRECISION A(LDA,*)
  62 C
  63       DOUBLE PRECISION T
  64       INTEGER IDAMAX,J,K,KP1,L,NM1
  65 C
  66 C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
  67 C
  68 C***FIRST EXECUTABLE STATEMENT  DGEFA
  69       INFO = 0
  70       NM1 = N - 1
  71       IF (NM1 .LT. 1) GO TO 70
  72       DO 60 K = 1, NM1
  73          KP1 = K + 1
  74 C
  75 C        FIND L = PIVOT INDEX
  76 C
  77          L = IDAMAX(N-K+1,A(K,K),1) + K - 1
  78          IPVT(K) = L
  79 C
  80 C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
  81 C
  82          IF (A(L,K) .EQ. 0.0D0) GO TO 40
  83 C
  84 C           INTERCHANGE IF NECESSARY
  85 C
  86             IF (L .EQ. K) GO TO 10
  87                T = A(L,K)
  88                A(L,K) = A(K,K)
  89                A(K,K) = T
  90    10       CONTINUE
  91 C
  92 C           COMPUTE MULTIPLIERS
  93 C
  94             T = -1.0D0/A(K,K)
  95             CALL DSCAL(N-K,T,A(K+1,K),1)
  96 C
  97 C           ROW ELIMINATION WITH COLUMN INDEXING
  98 C
  99             DO 30 J = KP1, N
 100                T = A(L,J)
 101                IF (L .EQ. K) GO TO 20
 102                   A(L,J) = A(K,J)
 103                   A(K,J) = T
 104    20          CONTINUE
 105                CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
 106    30       CONTINUE
 107          GO TO 50
 108    40    CONTINUE
 109             INFO = K
 110    50    CONTINUE
 111    60 CONTINUE
 112    70 CONTINUE
 113       IPVT(N) = N
 114       IF (A(N,N) .EQ. 0.0D0) INFO = N
 115       RETURN
 116       END