share/odepack/fortran/dgbsl.f

   1 *DECK DGBSL
   2       SUBROUTINE DGBSL (ABD, LDA, N, ML, MU, IPVT, B, JOB)
   3 C***BEGIN PROLOGUE  DGBSL
   4 C***PURPOSE  Solve the real band system A*X=B or TRANS(A)*X=B using
   5 C            the factors computed by DGBCO or DGBFA.
   6 C***CATEGORY  D2A2
   7 C***TYPE      DOUBLE PRECISION (SGBSL-S, DGBSL-D, CGBSL-C)
   8 C***KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
   9 C***AUTHOR  Moler, C. B., (U. of New Mexico)
  10 C***DESCRIPTION
  11 C
  12 C     DGBSL solves the double precision band system
  13 C     A * X = B  or  TRANS(A) * X = B
  14 C     using the factors computed by DGBCO or DGBFA.
  15 C
  16 C     On Entry
  17 C
  18 C        ABD     DOUBLE PRECISION(LDA, N)
  19 C                the output from DGBCO or DGBFA.
  20 C
  21 C        LDA     INTEGER
  22 C                the leading dimension of the array  ABD .
  23 C
  24 C        N       INTEGER
  25 C                the order of the original matrix.
  26 C
  27 C        ML      INTEGER
  28 C                number of diagonals below the main diagonal.
  29 C
  30 C        MU      INTEGER
  31 C                number of diagonals above the main diagonal.
  32 C
  33 C        IPVT    INTEGER(N)
  34 C                the pivot vector from DGBCO or DGBFA.
  35 C
  36 C        B       DOUBLE PRECISION(N)
  37 C                the right hand side vector.
  38 C
  39 C        JOB     INTEGER
  40 C                = 0         to solve  A*X = B ,
  41 C                = nonzero   to solve  TRANS(A)*X = B , where
  42 C                            TRANS(A)  is the transpose.
  43 C
  44 C     On Return
  45 C
  46 C        B       the solution vector  X .
  47 C
  48 C     Error Condition
  49 C
  50 C        A division by zero will occur if the input factor contains a
  51 C        zero on the diagonal.  Technically this indicates singularity
  52 C        but it is often caused by improper arguments or improper
  53 C        setting of LDA .  It will not occur if the subroutines are
  54 C        called correctly and if DGBCO has set RCOND .GT. 0.0
  55 C        or DGBFA has set INFO .EQ. 0 .
  56 C
  57 C     To compute  INVERSE(A) * C  where  C  is a matrix
  58 C     with  P  columns
  59 C           CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
  60 C           IF (RCOND is too small) GO TO ...
  61 C           DO 10 J = 1, P
  62 C              CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
  63 C        10 CONTINUE
  64 C
  65 C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
  66 C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
  67 C***ROUTINES CALLED  DAXPY, DDOT
  68 C***REVISION HISTORY  (YYMMDD)
  69 C   780814  DATE WRITTEN
  70 C   890531  Changed all specific intrinsics to generic.  (WRB)
  71 C   890831  Modified array declarations.  (WRB)
  72 C   890831  REVISION DATE from Version 3.2
  73 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  74 C   900326  Removed duplicate information from DESCRIPTION section.
  75 C           (WRB)
  76 C   920501  Reformatted the REFERENCES section.  (WRB)
  77 C***END PROLOGUE  DGBSL
  78       INTEGER LDA,N,ML,MU,IPVT(*),JOB
  79       DOUBLE PRECISION ABD(LDA,*),B(*)
  80 C
  81       DOUBLE PRECISION DDOT,T
  82       INTEGER K,KB,L,LA,LB,LM,M,NM1
  83 C***FIRST EXECUTABLE STATEMENT  DGBSL
  84       M = MU + ML + 1
  85       NM1 = N - 1
  86       IF (JOB .NE. 0) GO TO 50
  87 C
  88 C        JOB = 0 , SOLVE  A * X = B
  89 C        FIRST SOLVE L*Y = B
  90 C
  91          IF (ML .EQ. 0) GO TO 30
  92          IF (NM1 .LT. 1) GO TO 30
  93             DO 20 K = 1, NM1
  94                LM = MIN(ML,N-K)
  95                L = IPVT(K)
  96                T = B(L)
  97                IF (L .EQ. K) GO TO 10
  98                   B(L) = B(K)
  99                   B(K) = T
 100    10          CONTINUE
 101                CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
 102    20       CONTINUE
 103    30    CONTINUE
 104 C
 105 C        NOW SOLVE  U*X = Y
 106 C
 107          DO 40 KB = 1, N
 108             K = N + 1 - KB
 109             B(K) = B(K)/ABD(M,K)
 110             LM = MIN(K,M) - 1
 111             LA = M - LM
 112             LB = K - LM
 113             T = -B(K)
 114             CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
 115    40    CONTINUE
 116       GO TO 100
 117    50 CONTINUE
 118 C
 119 C        JOB = NONZERO, SOLVE  TRANS(A) * X = B
 120 C        FIRST SOLVE  TRANS(U)*Y = B
 121 C
 122          DO 60 K = 1, N
 123             LM = MIN(K,M) - 1
 124             LA = M - LM
 125             LB = K - LM
 126             T = DDOT(LM,ABD(LA,K),1,B(LB),1)
 127             B(K) = (B(K) - T)/ABD(M,K)
 128    60    CONTINUE
 129 C
 130 C        NOW SOLVE TRANS(L)*X = Y
 131 C
 132          IF (ML .EQ. 0) GO TO 90
 133          IF (NM1 .LT. 1) GO TO 90
 134             DO 80 KB = 1, NM1
 135                K = N - KB
 136                LM = MIN(ML,N-K)
 137                B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
 138                L = IPVT(K)
 139                IF (L .EQ. K) GO TO 70
 140                   T = B(L)
 141                   B(L) = B(K)
 142                   B(K) = T
 143    70          CONTINUE
 144    80       CONTINUE
 145    90    CONTINUE
 146   100 CONTINUE
 147       RETURN
 148       END