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