home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / linpack / sspsl.for < prev    next >
Encoding:
Text File  |  1984-01-07  |  4.4 KB  |  181 lines
  1.       SUBROUTINE SSPSL(AP,N,KPVT,B)
  2.       INTEGER N,KPVT(1)
  3.       REAL AP(1),B(1)
  4. C
  5. C     SSISL SOLVES THE REAL SYMMETRIC SYSTEM
  6. C     A * X = B
  7. C     USING THE FACTORS COMPUTED BY SSPFA.
  8. C
  9. C     ON ENTRY
  10. C
  11. C        AP      REAL(N*(N+1)/2)
  12. C                THE OUTPUT FROM SSPFA.
  13. C
  14. C        N       INTEGER
  15. C                THE ORDER OF THE MATRIX  A .
  16. C
  17. C        KPVT    INTEGER(N)
  18. C                THE PIVOT VECTOR FROM SSPFA.
  19. C
  20. C        B       REAL(N)
  21. C                THE RIGHT HAND SIDE VECTOR.
  22. C
  23. C     ON RETURN
  24. C
  25. C        B       THE SOLUTION VECTOR  X .
  26. C
  27. C     ERROR CONDITION
  28. C
  29. C        A DIVISION BY ZERO MAY OCCUR IF  SSPCO  HAS SET RCOND .EQ. 0.0
  30. C        OR  SSPFA  HAS SET INFO .NE. 0  .
  31. C
  32. C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
  33. C     WITH  P  COLUMNS
  34. C           CALL SSPFA(AP,N,KPVT,INFO)
  35. C           IF (INFO .NE. 0) GO TO ...
  36. C           DO 10 J = 1, P
  37. C              CALL SSPSL(AP,N,KPVT,C(1,J))
  38. C        10 CONTINUE
  39. C
  40. C     LINPACK. THIS VERSION DATED 08/14/78 .
  41. C     JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB.
  42. C
  43. C     SUBROUTINES AND FUNCTIONS
  44. C
  45. C     BLAS SAXPY,SDOT
  46. C     FORTRAN IABS
  47. C
  48. C     INTERNAL VARIABLES.
  49. C
  50.       REAL AK,AKM1,BK,BKM1,SDOT,DENOM,TEMP
  51.       INTEGER IK,IKM1,IKP1,K,KK,KM1K,KM1KM1,KP
  52. C
  53. C     LOOP BACKWARD APPLYING THE TRANSFORMATIONS AND
  54. C     D INVERSE TO B.
  55. C
  56.       K = N
  57.       IK = (N*(N - 1))/2
  58.    10 IF (K .EQ. 0) GO TO 80
  59.          KK = IK + K
  60.          IF (KPVT(K) .LT. 0) GO TO 40
  61. C
  62. C           1 X 1 PIVOT BLOCK.
  63. C
  64.             IF (K .EQ. 1) GO TO 30
  65.                KP = KPVT(K)
  66.                IF (KP .EQ. K) GO TO 20
  67. C
  68. C                 INTERCHANGE.
  69. C
  70.                   TEMP = B(K)
  71.                   B(K) = B(KP)
  72.                   B(KP) = TEMP
  73.    20          CONTINUE
  74. C
  75. C              APPLY THE TRANSFORMATION.
  76. C
  77.                CALL SAXPY(K-1,B(K),AP(IK+1),1,B(1),1)
  78.    30       CONTINUE
  79. C
  80. C           APPLY D INVERSE.
  81. C
  82.             B(K) = B(K)/AP(KK)
  83.             K = K - 1
  84.             IK = IK - K
  85.          GO TO 70
  86.    40    CONTINUE
  87. C
  88. C           2 X 2 PIVOT BLOCK.
  89. C
  90.             IKM1 = IK - (K - 1)
  91.             IF (K .EQ. 2) GO TO 60
  92.                KP = IABS(KPVT(K))
  93.                IF (KP .EQ. K - 1) GO TO 50
  94. C
  95. C                 INTERCHANGE.
  96. C
  97.                   TEMP = B(K-1)
  98.                   B(K-1) = B(KP)
  99.                   B(KP) = TEMP
  100.    50          CONTINUE
  101. C
  102. C              APPLY THE TRANSFORMATION.
  103. C
  104.                CALL SAXPY(K-2,B(K),AP(IK+1),1,B(1),1)
  105.                CALL SAXPY(K-2,B(K-1),AP(IKM1+1),1,B(1),1)
  106.    60       CONTINUE
  107. C
  108. C           APPLY D INVERSE.
  109. C
  110.             KM1K = IK + K - 1
  111.             KK = IK + K
  112.             AK = AP(KK)/AP(KM1K)
  113.             KM1KM1 = IKM1 + K - 1
  114.             AKM1 = AP(KM1KM1)/AP(KM1K)
  115.             BK = B(K)/AP(KM1K)
  116.             BKM1 = B(K-1)/AP(KM1K)
  117.             DENOM = AK*AKM1 - 1.0E0
  118.             B(K) = (AKM1*BK - BKM1)/DENOM
  119.             B(K-1) = (AK*BKM1 - BK)/DENOM
  120.             K = K - 2
  121.             IK = IK - (K + 1) - K
  122.    70    CONTINUE
  123.       GO TO 10
  124.    80 CONTINUE
  125. C
  126. C     LOOP FORWARD APPLYING THE TRANSFORMATIONS.
  127. C
  128.       K = 1
  129.       IK = 0
  130.    90 IF (K .GT. N) GO TO 160
  131.          IF (KPVT(K) .LT. 0) GO TO 120
  132. C
  133. C           1 X 1 PIVOT BLOCK.
  134. C
  135.             IF (K .EQ. 1) GO TO 110
  136. C
  137. C              APPLY THE TRANSFORMATION.
  138. C
  139.                B(K) = B(K) + SDOT(K-1,AP(IK+1),1,B(1),1)
  140.                KP = KPVT(K)
  141.                IF (KP .EQ. K) GO TO 100
  142. C
  143. C                 INTERCHANGE.
  144. C
  145.                   TEMP = B(K)
  146.                   B(K) = B(KP)
  147.                   B(KP) = TEMP
  148.   100          CONTINUE
  149.   110       CONTINUE
  150.             IK = IK + K
  151.             K = K + 1
  152.          GO TO 150
  153.   120    CONTINUE
  154. C
  155. C           2 X 2 PIVOT BLOCK.
  156. C
  157.             IF (K .EQ. 1) GO TO 140
  158. C
  159. C              APPLY THE TRANSFORMATION.
  160. C
  161.                B(K) = B(K) + SDOT(K-1,AP(IK+1),1,B(1),1)
  162.                IKP1 = IK + K
  163.                B(K+1) = B(K+1) + SDOT(K-1,AP(IKP1+1),1,B(1),1)
  164.                KP = IABS(KPVT(K))
  165.                IF (KP .EQ. K) GO TO 130
  166. C
  167. C                 INTERCHANGE.
  168. C
  169.                   TEMP = B(K)
  170.                   B(K) = B(KP)
  171.                   B(KP) = TEMP
  172.   130          CONTINUE
  173.   140       CONTINUE
  174.             IK = IK + K + K + 1
  175.             K = K + 2
  176.   150    CONTINUE
  177.       GO TO 90
  178.   160 CONTINUE
  179.       RETURN
  180.       END
  181.