home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / fortran / linpklib.arc / SSPFA.FOR < prev    next >
Encoding:
Text File  |  1984-01-08  |  7.8 KB  |  265 lines
  1.       SUBROUTINE SSPFA(AP,N,KPVT,INFO)
  2.       INTEGER N,KPVT(1),INFO
  3.       REAL AP(1)
  4. C
  5. C     SSPFA FACTORS A REAL SYMMETRIC MATRIX STORED IN
  6. C     PACKED FORM BY ELIMINATION WITH SYMMETRIC PIVOTING.
  7. C
  8. C     TO SOLVE  A*X = B , FOLLOW SSPFA BY SSPSL.
  9. C     TO COMPUTE  INVERSE(A)*C , FOLLOW SSPFA BY SSPSL.
  10. C     TO COMPUTE  DETERMINANT(A) , FOLLOW SSPFA BY SSPDI.
  11. C     TO COMPUTE  INERTIA(A) , FOLLOW SSPFA BY SSPDI.
  12. C     TO COMPUTE  INVERSE(A) , FOLLOW SSPFA BY SSPDI.
  13. C
  14. C     ON ENTRY
  15. C
  16. C        AP      REAL (N*(N+1)/2)
  17. C                THE PACKED FORM OF A SYMMETRIC MATRIX  A .  THE
  18. C                COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY
  19. C                IN A ONE-DIMENSIONAL ARRAY OF LENGTH  N*(N+1)/2 .
  20. C                SEE COMMENTS BELOW FOR DETAILS.
  21. C
  22. C        N       INTEGER
  23. C                THE ORDER OF THE MATRIX  A .
  24. C
  25. C     OUTPUT
  26. C
  27. C        AP      A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH
  28. C                WERE USED TO OBTAIN IT STORED IN PACKED FORM.
  29. C                THE FACTORIZATION CAN BE WRITTEN  A = U*D*TRANS(U)
  30. C                WHERE  U  IS A PRODUCT OF PERMUTATION AND UNIT
  31. C                UPPER TRIANGULAR MATRICES , TRANS(U) IS THE
  32. C                TRANSPOSE OF  U , AND  D  IS BLOCK DIAGONAL
  33. C                WITH 1 BY 1 AND 2 BY 2 BLOCKS.
  34. C
  35. C        KPVT    INTEGER(N)
  36. C                AN INTEGER VECTOR OF PIVOT INDICES.
  37. C
  38. C        INFO    INTEGER
  39. C                = 0  NORMAL VALUE.
  40. C                = K  IF THE K-TH PIVOT BLOCK IS SINGULAR. THIS IS
  41. C                     NOT AN ERROR CONDITION FOR THIS SUBROUTINE,
  42. C                     BUT IT DOES INDICATE THAT SSPSL OR SSPDI MAY
  43. C                     DIVIDE BY ZERO IF CALLED.
  44. C
  45. C     PACKED STORAGE
  46. C
  47. C          THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER
  48. C          TRIANGLE OF A SYMMETRIC MATRIX.
  49. C
  50. C                K = 0
  51. C                DO 20 J = 1, N
  52. C                   DO 10 I = 1, J
  53. C                      K = K + 1
  54. C                      AP(K)  = A(I,J)
  55. C             10    CONTINUE
  56. C             20 CONTINUE
  57. C
  58. C     LINPACK. THIS VERSION DATED 08/14/78 .
  59. C     JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB.
  60. C
  61. C     SUBROUTINES AND FUNCTIONS
  62. C
  63. C     BLAS SAXPY,SSWAP,ISAMAX
  64. C     FORTRAN ABS,AMAX1,SQRT
  65. C
  66. C     INTERNAL VARIABLES
  67. C
  68.       REAL AK,AKM1,BK,BKM1,DENOM,MULK,MULKM1,T
  69.       REAL ABSAKK,ALPHA,COLMAX,ROWMAX
  70.       INTEGER ISAMAX,IJ,IJJ,IK,IKM1,IM,IMAX,IMAXP1,IMIM,IMJ,IMK
  71.       INTEGER J,JJ,JK,JKM1,JMAX,JMIM,K,KK,KM1,KM1K,KM1KM1,KM2,KSTEP
  72.       LOGICAL SWAP
  73. C
  74. C
  75. C     INITIALIZE
  76. C
  77. C     ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE.
  78.       ALPHA = (1.0E0 + SQRT(17.0E0))/8.0E0
  79. C
  80.       INFO = 0
  81. C
  82. C     MAIN LOOP ON K, WHICH GOES FROM N TO 1.
  83. C
  84.       K = N
  85.       IK = (N*(N - 1))/2
  86.    10 CONTINUE
  87. C
  88. C        LEAVE THE LOOP IF K=0 OR K=1.
  89. C
  90. C     ...EXIT
  91.          IF (K .EQ. 0) GO TO 200
  92.          IF (K .GT. 1) GO TO 20
  93.             KPVT(1) = 1
  94.             IF (AP(1) .EQ. 0.0E0) INFO = 1
  95. C     ......EXIT
  96.             GO TO 200
  97.    20    CONTINUE
  98. C
  99. C        THIS SECTION OF CODE DETERMINES THE KIND OF
  100. C        ELIMINATION TO BE PERFORMED.  WHEN IT IS COMPLETED,
  101. C        KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND
  102. C        SWAP WILL BE SET TO .TRUE. IF AN INTERCHANGE IS
  103. C        REQUIRED.
  104. C
  105.          KM1 = K - 1
  106.          KK = IK + K
  107.          ABSAKK = ABS(AP(KK))
  108. C
  109. C        DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
  110. C        COLUMN K.
  111. C
  112.          IMAX = ISAMAX(K-1,AP(IK+1),1)
  113.          IMK = IK + IMAX
  114.          COLMAX = ABS(AP(IMK))
  115.          IF (ABSAKK .LT. ALPHA*COLMAX) GO TO 30
  116.             KSTEP = 1
  117.             SWAP = .FALSE.
  118.          GO TO 90
  119.    30    CONTINUE
  120. C
  121. C           DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
  122. C           ROW IMAX.
  123. C
  124.             ROWMAX = 0.0E0
  125.             IMAXP1 = IMAX + 1
  126.             IM = IMAX*(IMAX - 1)/2
  127.             IMJ = IM + 2*IMAX
  128.             DO 40 J = IMAXP1, K
  129.                ROWMAX = AMAX1(ROWMAX,ABS(AP(IMJ)))
  130.                IMJ = IMJ + J
  131.    40       CONTINUE
  132.             IF (IMAX .EQ. 1) GO TO 50
  133.                JMAX = ISAMAX(IMAX-1,AP(IM+1),1)
  134.                JMIM = JMAX + IM
  135.                ROWMAX = AMAX1(ROWMAX,ABS(AP(JMIM)))
  136.    50       CONTINUE
  137.             IMIM = IMAX + IM
  138.             IF (ABS(AP(IMIM)) .LT. ALPHA*ROWMAX) GO TO 60
  139.                KSTEP = 1
  140.                SWAP = .TRUE.
  141.             GO TO 80
  142.    60       CONTINUE
  143.             IF (ABSAKK .LT. ALPHA*COLMAX*(COLMAX/ROWMAX)) GO TO 70
  144.                KSTEP = 1
  145.                SWAP = .FALSE.
  146.             GO TO 80
  147.    70       CONTINUE
  148.                KSTEP = 2
  149.                SWAP = IMAX .NE. KM1
  150.    80       CONTINUE
  151.    90    CONTINUE
  152.          IF (AMAX1(ABSAKK,COLMAX) .NE. 0.0E0) GO TO 100
  153. C
  154. C           COLUMN K IS ZERO.  SET INFO AND ITERATE THE LOOP.
  155. C
  156.             KPVT(K) = K
  157.             INFO = K
  158.          GO TO 190
  159.   100    CONTINUE
  160.          IF (KSTEP .EQ. 2) GO TO 140
  161. C
  162. C           1 X 1 PIVOT BLOCK.
  163. C
  164.             IF (.NOT.SWAP) GO TO 120
  165. C
  166. C              PERFORM AN INTERCHANGE.
  167. C
  168.                CALL SSWAP(IMAX,AP(IM+1),1,AP(IK+1),1)
  169.                IMJ = IK + IMAX
  170.                DO 110 JJ = IMAX, K
  171.                   J = K + IMAX - JJ
  172.                   JK = IK + J
  173.                   T = AP(JK)
  174.                   AP(JK) = AP(IMJ)
  175.                   AP(IMJ) = T
  176.                   IMJ = IMJ - (J - 1)
  177.   110          CONTINUE
  178.   120       CONTINUE
  179. C
  180. C           PERFORM THE ELIMINATION.
  181. C
  182.             IJ = IK - (K - 1)
  183.             DO 130 JJ = 1, KM1
  184.                J = K - JJ
  185.                JK = IK + J
  186.                MULK = -AP(JK)/AP(KK)
  187.                T = MULK
  188.                CALL SAXPY(J,T,AP(IK+1),1,AP(IJ+1),1)
  189.                IJJ = IJ + J
  190.                AP(JK) = MULK
  191.                IJ = IJ - (J - 1)
  192.   130       CONTINUE
  193. C
  194. C           SET THE PIVOT ARRAY.
  195. C
  196.             KPVT(K) = K
  197.             IF (SWAP) KPVT(K) = IMAX
  198.          GO TO 190
  199.   140    CONTINUE
  200. C
  201. C           2 X 2 PIVOT BLOCK.
  202. C
  203.             KM1K = IK + K - 1
  204.             IKM1 = IK - (K - 1)
  205.             IF (.NOT.SWAP) GO TO 160
  206. C
  207. C              PERFORM AN INTERCHANGE.
  208. C
  209.                CALL SSWAP(IMAX,AP(IM+1),1,AP(IKM1+1),1)
  210.                IMJ = IKM1 + IMAX
  211.                DO 150 JJ = IMAX, KM1
  212.                   J = KM1 + IMAX - JJ
  213.                   JKM1 = IKM1 + J
  214.                   T = AP(JKM1)
  215.                   AP(JKM1) = AP(IMJ)
  216.                   AP(IMJ) = T
  217.                   IMJ = IMJ - (J - 1)
  218.   150          CONTINUE
  219.                T = AP(KM1K)
  220.                AP(KM1K) = AP(IMK)
  221.                AP(IMK) = T
  222.   160       CONTINUE
  223. C
  224. C           PERFORM THE ELIMINATION.
  225. C
  226.             KM2 = K - 2
  227.             IF (KM2 .EQ. 0) GO TO 180
  228.                AK = AP(KK)/AP(KM1K)
  229.                KM1KM1 = IKM1 + K - 1
  230.                AKM1 = AP(KM1KM1)/AP(KM1K)
  231.                DENOM = 1.0E0 - AK*AKM1
  232.                IJ = IK - (K - 1) - (K - 2)
  233.                DO 170 JJ = 1, KM2
  234.                   J = KM1 - JJ
  235.                   JK = IK + J
  236.                   BK = AP(JK)/AP(KM1K)
  237.                   JKM1 = IKM1 + J
  238.                   BKM1 = AP(JKM1)/AP(KM1K)
  239.                   MULK = (AKM1*BK - BKM1)/DENOM
  240.                   MULKM1 = (AK*BKM1 - BK)/DENOM
  241.                   T = MULK
  242.                   CALL SAXPY(J,T,AP(IK+1),1,AP(IJ+1),1)
  243.                   T = MULKM1
  244.                   CALL SAXPY(J,T,AP(IKM1+1),1,AP(IJ+1),1)
  245.                   AP(JK) = MULK
  246.                   AP(JKM1) = MULKM1
  247.                   IJJ = IJ + J
  248.                   IJ = IJ - (J - 1)
  249.   170          CONTINUE
  250.   180       CONTINUE
  251. C
  252. C           SET THE PIVOT ARRAY.
  253. C
  254.             KPVT(K) = 1 - K
  255.             IF (SWAP) KPVT(K) = -IMAX
  256.             KPVT(K-1) = KPVT(K)
  257.   190    CONTINUE
  258.          IK = IK - (K - 1)
  259.          IF (KSTEP .EQ. 2) IK = IK - (K - 2)
  260.          K = K - KSTEP
  261.       GO TO 10
  262.   200 CONTINUE
  263.       RETURN
  264.       END
  265.