share/odepack/fortran/dgbfa.f

   1 *DECK DGBFA
   2       SUBROUTINE DGBFA (ABD, LDA, N, ML, MU, IPVT, INFO)
   3 C***BEGIN PROLOGUE  DGBFA
   4 C***PURPOSE  Factor a band matrix using Gaussian elimination.
   5 C***CATEGORY  D2A2
   6 C***TYPE      DOUBLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
   7 C***KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
   8 C***AUTHOR  Moler, C. B., (U. of New Mexico)
   9 C***DESCRIPTION
  10 C
  11 C     DGBFA factors a double precision band matrix by elimination.
  12 C
  13 C     DGBFA is usually called by DGBCO, but it can be called
  14 C     directly with a saving in time if  RCOND  is not needed.
  15 C
  16 C     On Entry
  17 C
  18 C        ABD     DOUBLE PRECISION(LDA, N)
  19 C                contains the matrix in band storage.  The columns
  20 C                of the matrix are stored in the columns of  ABD  and
  21 C                the diagonals of the matrix are stored in rows
  22 C                ML+1 through 2*ML+MU+1 of  ABD .
  23 C                See the comments below for details.
  24 C
  25 C        LDA     INTEGER
  26 C                the leading dimension of the array  ABD .
  27 C                LDA must be .GE. 2*ML + MU + 1 .
  28 C
  29 C        N       INTEGER
  30 C                the order of the original matrix.
  31 C
  32 C        ML      INTEGER
  33 C                number of diagonals below the main diagonal.
  34 C                0 .LE. ML .LT.  N .
  35 C
  36 C        MU      INTEGER
  37 C                number of diagonals above the main diagonal.
  38 C                0 .LE. MU .LT.  N .
  39 C                More efficient if  ML .LE. MU .
  40 C     On Return
  41 C
  42 C        ABD     an upper triangular matrix in band storage and
  43 C                the multipliers which were used to obtain it.
  44 C                The factorization can be written  A = L*U  where
  45 C                L  is a product of permutation and unit lower
  46 C                triangular matrices and  U  is upper triangular.
  47 C
  48 C        IPVT    INTEGER(N)
  49 C                an integer vector of pivot indices.
  50 C
  51 C        INFO    INTEGER
  52 C                = 0  normal value.
  53 C                = K  if  U(K,K) .EQ. 0.0 .  This is not an error
  54 C                     condition for this subroutine, but it does
  55 C                     indicate that DGBSL will divide by zero if
  56 C                     called.  Use  RCOND  in DGBCO for a reliable
  57 C                     indication of singularity.
  58 C
  59 C     Band Storage
  60 C
  61 C           If  A  is a band matrix, the following program segment
  62 C           will set up the input.
  63 C
  64 C                   ML = (band width below the diagonal)
  65 C                   MU = (band width above the diagonal)
  66 C                   M = ML + MU + 1
  67 C                   DO 20 J = 1, N
  68 C                      I1 = MAX(1, J-MU)
  69 C                      I2 = MIN(N, J+ML)
  70 C                      DO 10 I = I1, I2
  71 C                         K = I - J + M
  72 C                         ABD(K,J) = A(I,J)
  73 C                10    CONTINUE
  74 C                20 CONTINUE
  75 C
  76 C           This uses rows  ML+1  through  2*ML+MU+1  of  ABD .
  77 C           In addition, the first  ML  rows in  ABD  are used for
  78 C           elements generated during the triangularization.
  79 C           The total number of rows needed in  ABD  is  2*ML+MU+1 .
  80 C           The  ML+MU by ML+MU  upper left triangle and the
  81 C           ML by ML  lower right triangle are not referenced.
  82 C
  83 C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
  84 C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
  85 C***ROUTINES CALLED  DAXPY, DSCAL, IDAMAX
  86 C***REVISION HISTORY  (YYMMDD)
  87 C   780814  DATE WRITTEN
  88 C   890531  Changed all specific intrinsics to generic.  (WRB)
  89 C   890831  Modified array declarations.  (WRB)
  90 C   890831  REVISION DATE from Version 3.2
  91 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  92 C   900326  Removed duplicate information from DESCRIPTION section.
  93 C           (WRB)
  94 C   920501  Reformatted the REFERENCES section.  (WRB)
  95 C***END PROLOGUE  DGBFA
  96       INTEGER LDA,N,ML,MU,IPVT(*),INFO
  97       DOUBLE PRECISION ABD(LDA,*)
  98 C
  99       DOUBLE PRECISION T
 100       INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
 101 C
 102 C***FIRST EXECUTABLE STATEMENT  DGBFA
 103       M = ML + MU + 1
 104       INFO = 0
 105 C
 106 C     ZERO INITIAL FILL-IN COLUMNS
 107 C
 108       J0 = MU + 2
 109       J1 = MIN(N,M) - 1
 110       IF (J1 .LT. J0) GO TO 30
 111       DO 20 JZ = J0, J1
 112          I0 = M + 1 - JZ
 113          DO 10 I = I0, ML
 114             ABD(I,JZ) = 0.0D0
 115    10    CONTINUE
 116    20 CONTINUE
 117    30 CONTINUE
 118       JZ = J1
 119       JU = 0
 120 C
 121 C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
 122 C
 123       NM1 = N - 1
 124       IF (NM1 .LT. 1) GO TO 130
 125       DO 120 K = 1, NM1
 126          KP1 = K + 1
 127 C
 128 C        ZERO NEXT FILL-IN COLUMN
 129 C
 130          JZ = JZ + 1
 131          IF (JZ .GT. N) GO TO 50
 132          IF (ML .LT. 1) GO TO 50
 133             DO 40 I = 1, ML
 134                ABD(I,JZ) = 0.0D0
 135    40       CONTINUE
 136    50    CONTINUE
 137 C
 138 C        FIND L = PIVOT INDEX
 139 C
 140          LM = MIN(ML,N-K)
 141          L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
 142          IPVT(K) = L + K - M
 143 C
 144 C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
 145 C
 146          IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
 147 C
 148 C           INTERCHANGE IF NECESSARY
 149 C
 150             IF (L .EQ. M) GO TO 60
 151                T = ABD(L,K)
 152                ABD(L,K) = ABD(M,K)
 153                ABD(M,K) = T
 154    60       CONTINUE
 155 C
 156 C           COMPUTE MULTIPLIERS
 157 C
 158             T = -1.0D0/ABD(M,K)
 159             CALL DSCAL(LM,T,ABD(M+1,K),1)
 160 C
 161 C           ROW ELIMINATION WITH COLUMN INDEXING
 162 C
 163             JU = MIN(MAX(JU,MU+IPVT(K)),N)
 164             MM = M
 165             IF (JU .LT. KP1) GO TO 90
 166             DO 80 J = KP1, JU
 167                L = L - 1
 168                MM = MM - 1
 169                T = ABD(L,J)
 170                IF (L .EQ. MM) GO TO 70
 171                   ABD(L,J) = ABD(MM,J)
 172                   ABD(MM,J) = T
 173    70          CONTINUE
 174                CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
 175    80       CONTINUE
 176    90       CONTINUE
 177          GO TO 110
 178   100    CONTINUE
 179             INFO = K
 180   110    CONTINUE
 181   120 CONTINUE
 182   130 CONTINUE
 183       IPVT(N) = N
 184       IF (ABD(M,N) .EQ. 0.0D0) INFO = N
 185       RETURN
 186       END