home *** CD-ROM | disk | FTP | other *** search
- SUBROUTINE SSISL(A,LDA,N,KPVT,B)
- INTEGER LDA,N,KPVT(1)
- REAL A(LDA,1),B(1)
- C
- C SSISL SOLVES THE REAL SYMMETRIC SYSTEM
- C A * X = B
- C USING THE FACTORS COMPUTED BY SSIFA.
- C
- C ON ENTRY
- C
- C A REAL(LDA,N)
- C THE OUTPUT FROM SSIFA.
- C
- C LDA INTEGER
- C THE LEADING DIMENSION OF THE ARRAY A .
- C
- C N INTEGER
- C THE ORDER OF THE MATRIX A .
- C
- C KPVT INTEGER(N)
- C THE PIVOT VECTOR FROM SSIFA.
- C
- C B REAL(N)
- C THE RIGHT HAND SIDE VECTOR.
- C
- C ON RETURN
- C
- C B THE SOLUTION VECTOR X .
- C
- C ERROR CONDITION
- C
- C A DIVISION BY ZERO MAY OCCUR IF SSICO HAS SET RCOND .EQ. 0.0
- C OR SSIFA HAS SET INFO .NE. 0 .
- C
- C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
- C WITH P COLUMNS
- C CALL SSIFA(A,LDA,N,KPVT,INFO)
- C IF (INFO .NE. 0) GO TO ...
- C DO 10 J = 1, P
- C CALL SSISL(A,LDA,N,KPVT,C(1,J))
- C 10 CONTINUE
- C
- C LINPACK. THIS VERSION DATED 08/14/78 .
- C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB.
- C
- C SUBROUTINES AND FUNCTIONS
- C
- C BLAS SAXPY,SDOT
- C FORTRAN IABS
- C
- C INTERNAL VARIABLES.
- C
- REAL AK,AKM1,BK,BKM1,SDOT,DENOM,TEMP
- INTEGER K,KP
- C
- C LOOP BACKWARD APPLYING THE TRANSFORMATIONS AND
- C D INVERSE TO B.
- C
- K = N
- 10 IF (K .EQ. 0) GO TO 80
- IF (KPVT(K) .LT. 0) GO TO 40
- C
- C 1 X 1 PIVOT BLOCK.
- C
- IF (K .EQ. 1) GO TO 30
- KP = KPVT(K)
- IF (KP .EQ. K) GO TO 20
- C
- C INTERCHANGE.
- C
- TEMP = B(K)
- B(K) = B(KP)
- B(KP) = TEMP
- 20 CONTINUE
- C
- C APPLY THE TRANSFORMATION.
- C
- CALL SAXPY(K-1,B(K),A(1,K),1,B(1),1)
- 30 CONTINUE
- C
- C APPLY D INVERSE.
- C
- B(K) = B(K)/A(K,K)
- K = K - 1
- GO TO 70
- 40 CONTINUE
- C
- C 2 X 2 PIVOT BLOCK.
- C
- IF (K .EQ. 2) GO TO 60
- KP = IABS(KPVT(K))
- IF (KP .EQ. K - 1) GO TO 50
- C
- C INTERCHANGE.
- C
- TEMP = B(K-1)
- B(K-1) = B(KP)
- B(KP) = TEMP
- 50 CONTINUE
- C
- C APPLY THE TRANSFORMATION.
- C
- CALL SAXPY(K-2,B(K),A(1,K),1,B(1),1)
- CALL SAXPY(K-2,B(K-1),A(1,K-1),1,B(1),1)
- 60 CONTINUE
- C
- C APPLY D INVERSE.
- C
- AK = A(K,K)/A(K-1,K)
- AKM1 = A(K-1,K-1)/A(K-1,K)
- BK = B(K)/A(K-1,K)
- BKM1 = B(K-1)/A(K-1,K)
- DENOM = AK*AKM1 - 1.0E0
- B(K) = (AKM1*BK - BKM1)/DENOM
- B(K-1) = (AK*BKM1 - BK)/DENOM
- K = K - 2
- 70 CONTINUE
- GO TO 10
- 80 CONTINUE
- C
- C LOOP FORWARD APPLYING THE TRANSFORMATIONS.
- C
- K = 1
- 90 IF (K .GT. N) GO TO 160
- IF (KPVT(K) .LT. 0) GO TO 120
- C
- C 1 X 1 PIVOT BLOCK.
- C
- IF (K .EQ. 1) GO TO 110
- C
- C APPLY THE TRANSFORMATION.
- C
- B(K) = B(K) + SDOT(K-1,A(1,K),1,B(1),1)
- KP = KPVT(K)
- IF (KP .EQ. K) GO TO 100
- C
- C INTERCHANGE.
- C
- TEMP = B(K)
- B(K) = B(KP)
- B(KP) = TEMP
- 100 CONTINUE
- 110 CONTINUE
- K = K + 1
- GO TO 150
- 120 CONTINUE
- C
- C 2 X 2 PIVOT BLOCK.
- C
- IF (K .EQ. 1) GO TO 140
- C
- C APPLY THE TRANSFORMATION.
- C
- B(K) = B(K) + SDOT(K-1,A(1,K),1,B(1),1)
- B(K+1) = B(K+1) + SDOT(K-1,A(1,K+1),1,B(1),1)
- KP = IABS(KPVT(K))
- IF (KP .EQ. K) GO TO 130
- C
- C INTERCHANGE.
- C
- TEMP = B(K)
- B(K) = B(KP)
- B(KP) = TEMP
- 130 CONTINUE
- 140 CONTINUE
- K = K + 2
- 150 CONTINUE
- GO TO 90
- 160 CONTINUE
- RETURN
- END