| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149 |
- *DECK SGBSL
- SUBROUTINE SGBSL (ABD, LDA, N, ML, MU, IPVT, B, JOB)
- C***BEGIN PROLOGUE SGBSL
- C***PURPOSE Solve the real band system A*X=B or TRANS(A)*X=B using
- C the factors computed by SGBCO or SGBFA.
- C***LIBRARY SLATEC (LINPACK)
- C***CATEGORY D2A2
- C***TYPE SINGLE PRECISION (SGBSL-S, DGBSL-D, CGBSL-C)
- C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
- C***AUTHOR Moler, C. B., (U. of New Mexico)
- C***DESCRIPTION
- C
- C SGBSL solves the real band system
- C A * X = B or TRANS(A) * X = B
- C using the factors computed by SBGCO or SGBFA.
- C
- C On Entry
- C
- C ABD REAL(LDA, N)
- C the output from SBGCO or SGBFA.
- C
- C LDA INTEGER
- C the leading dimension of the array ABD .
- C
- C N INTEGER
- C the order of the original matrix.
- C
- C ML INTEGER
- C number of diagonals below the main diagonal.
- C
- C MU INTEGER
- C number of diagonals above the main diagonal.
- C
- C IPVT INTEGER(N)
- C the pivot vector from SBGCO or SGBFA.
- C
- C B REAL(N)
- C the right hand side vector.
- C
- C JOB INTEGER
- C = 0 to solve A*X = B ,
- C = nonzero to solve TRANS(A)*X = B , where
- C TRANS(A) is the transpose.
- C
- C On Return
- C
- C B the solution vector X .
- C
- C Error Condition
- C
- C A division by zero will occur if the input factor contains a
- C zero on the diagonal. Technically, this indicates singularity,
- C but it is often caused by improper arguments or improper
- C setting of LDA . It will not occur if the subroutines are
- C called correctly and if SBGCO has set RCOND .GT. 0.0
- C or SGBFA has set INFO .EQ. 0 .
- C
- C To compute INVERSE(A) * C where C is a matrix
- C with P columns
- C CALL SBGCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
- C If (RCOND is too small) GO TO ...
- C DO 10 J = 1, P
- C CALL SGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
- C 10 CONTINUE
- 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 SAXPY, SDOT
- C***REVISION HISTORY (YYMMDD)
- C 780814 DATE WRITTEN
- C 890531 Changed all specific intrinsics to generic. (WRB)
- 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 SGBSL
- INTEGER LDA,N,ML,MU,IPVT(*),JOB
- REAL ABD(LDA,*),B(*)
- C
- REAL SDOT,T
- INTEGER K,KB,L,LA,LB,LM,M,NM1
- C***FIRST EXECUTABLE STATEMENT SGBSL
- M = MU + ML + 1
- NM1 = N - 1
- IF (JOB .NE. 0) GO TO 50
- C
- C JOB = 0 , SOLVE A * X = B
- C FIRST SOLVE L*Y = B
- C
- IF (ML .EQ. 0) GO TO 30
- IF (NM1 .LT. 1) GO TO 30
- DO 20 K = 1, NM1
- LM = MIN(ML,N-K)
- L = IPVT(K)
- T = B(L)
- IF (L .EQ. K) GO TO 10
- B(L) = B(K)
- B(K) = T
- 10 CONTINUE
- CALL SAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
- 20 CONTINUE
- 30 CONTINUE
- C
- C NOW SOLVE U*X = Y
- C
- DO 40 KB = 1, N
- K = N + 1 - KB
- B(K) = B(K)/ABD(M,K)
- LM = MIN(K,M) - 1
- LA = M - LM
- LB = K - LM
- T = -B(K)
- CALL SAXPY(LM,T,ABD(LA,K),1,B(LB),1)
- 40 CONTINUE
- GO TO 100
- 50 CONTINUE
- C
- C JOB = NONZERO, SOLVE TRANS(A) * X = B
- C FIRST SOLVE TRANS(U)*Y = B
- C
- DO 60 K = 1, N
- LM = MIN(K,M) - 1
- LA = M - LM
- LB = K - LM
- T = SDOT(LM,ABD(LA,K),1,B(LB),1)
- B(K) = (B(K) - T)/ABD(M,K)
- 60 CONTINUE
- C
- C NOW SOLVE TRANS(L)*X = Y
- C
- IF (ML .EQ. 0) GO TO 90
- IF (NM1 .LT. 1) GO TO 90
- DO 80 KB = 1, NM1
- K = N - KB
- LM = MIN(ML,N-K)
- B(K) = B(K) + SDOT(LM,ABD(M+1,K),1,B(K+1),1)
- L = IPVT(K)
- IF (L .EQ. K) GO TO 70
- T = B(L)
- B(L) = B(K)
- B(K) = T
- 70 CONTINUE
- 80 CONTINUE
- 90 CONTINUE
- 100 CONTINUE
- RETURN
- END
|
