share/odepack/fortran/dgesl.f

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