home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / linpack / sspco.for < prev    next >
Encoding:
Text File  |  1984-01-07  |  8.7 KB  |  286 lines
  1.       SUBROUTINE SSPCO(AP,N,KPVT,RCOND,Z)
  2.       INTEGER N,KPVT(1)
  3.       REAL AP(1),Z(1)
  4.       REAL RCOND
  5. C
  6. C     SSPCO FACTORS A REAL SYMMETRIC MATRIX STORED IN PACKED
  7. C     FORM BY ELIMINATION WITH SYMMETRIC PIVOTING AND ESTIMATES
  8. C     THE CONDITION OF THE MATRIX.
  9. C
  10. C     IF  RCOND  IS NOT NEEDED, SSPFA IS SLIGHTLY FASTER.
  11. C     TO SOLVE  A*X = B , FOLLOW SSPCO BY SSPSL.
  12. C     TO COMPUTE  INVERSE(A)*C , FOLLOW SSPCO BY SSPSL.
  13. C     TO COMPUTE  INVERSE(A) , FOLLOW SSPCO BY SSPDI.
  14. C     TO COMPUTE  DETERMINANT(A) , FOLLOW SSPCO BY SSPDI.
  15. C     TO COMPUTE  INERTIA(A), FOLLOW SSPCO BY SSPDI.
  16. C
  17. C     ON ENTRY
  18. C
  19. C        AP      REAL (N*(N+1)/2)
  20. C                THE PACKED FORM OF A SYMMETRIC MATRIX  A .  THE
  21. C                COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY
  22. C                IN A ONE-DIMENSIONAL ARRAY OF LENGTH  N*(N+1)/2 .
  23. C                SEE COMMENTS BELOW FOR DETAILS.
  24. C
  25. C        N       INTEGER
  26. C                THE ORDER OF THE MATRIX  A .
  27. C
  28. C     OUTPUT
  29. C
  30. C        AP      A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH
  31. C                WERE USED TO OBTAIN IT STORED IN PACKED FORM.
  32. C                THE FACTORIZATION CAN BE WRITTEN  A = U*D*TRANS(U)
  33. C                WHERE  U  IS A PRODUCT OF PERMUTATION AND UNIT
  34. C                UPPER TRIANGULAR MATRICES , TRANS(U) IS THE
  35. C                TRANSPOSE OF  U , AND  D  IS BLOCK DIAGONAL
  36. C                WITH 1 BY 1 AND 2 BY 2 BLOCKS.
  37. C
  38. C        KPVT    INTEGER(N)
  39. C                AN INTEGER VECTOR OF PIVOT INDICES.
  40. C
  41. C        RCOND   REAL
  42. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .
  43. C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS
  44. C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE
  45. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
  46. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
  47. C                           1.0 + RCOND .EQ. 1.0
  48. C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING
  49. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
  50. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
  51. C                UNDERFLOWS.
  52. C
  53. C        Z       REAL(N)
  54. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
  55. C                IF  A  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS
  56. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
  57. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
  58. C
  59. C     PACKED STORAGE
  60. C
  61. C          THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER
  62. C          TRIANGLE OF A SYMMETRIC MATRIX.
  63. C
  64. C                K = 0
  65. C                DO 20 J = 1, N
  66. C                   DO 10 I = 1, J
  67. C                      K = K + 1
  68. C                      AP(K) = A(I,J)
  69. C             10    CONTINUE
  70. C             20 CONTINUE
  71. C
  72. C     LINPACK. THIS VERSION DATED 08/14/78 .
  73. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  74. C
  75. C     SUBROUTINES AND FUNCTIONS
  76. C
  77. C     LINPACK SSPFA
  78. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  79. C     FORTRAN ABS,AMAX1,IABS,SIGN
  80. C
  81. C     INTERNAL VARIABLES
  82. C
  83.       REAL AK,AKM1,BK,BKM1,SDOT,DENOM,EK,T
  84.       REAL ANORM,S,SASUM,YNORM
  85.       INTEGER I,IJ,IK,IKM1,IKP1,INFO,J,JM1,J1
  86.       INTEGER K,KK,KM1K,KM1KM1,KP,KPS,KS
  87. C
  88. C
  89. C     FIND NORM OF A USING ONLY UPPER HALF
  90. C
  91.       J1 = 1
  92.       DO 30 J = 1, N
  93.          Z(J) = SASUM(J,AP(J1),1)
  94.          IJ = J1
  95.          J1 = J1 + J
  96.          JM1 = J - 1
  97.          IF (JM1 .LT. 1) GO TO 20
  98.          DO 10 I = 1, JM1
  99.             Z(I) = Z(I) + ABS(AP(IJ))
  100.             IJ = IJ + 1
  101.    10    CONTINUE
  102.    20    CONTINUE
  103.    30 CONTINUE
  104.       ANORM = 0.0E0
  105.       DO 40 J = 1, N
  106.          ANORM = AMAX1(ANORM,Z(J))
  107.    40 CONTINUE
  108. C
  109. C     FACTOR
  110. C
  111.       CALL SSPFA(AP,N,KPVT,INFO)
  112. C
  113. C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
  114. C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
  115. C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
  116. C     GROWTH IN THE ELEMENTS OF W  WHERE  U*D*W = E .
  117. C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
  118. C
  119. C     SOLVE U*D*W = E
  120. C
  121.       EK = 1.0E0
  122.       DO 50 J = 1, N
  123.          Z(J) = 0.0E0
  124.    50 CONTINUE
  125.       K = N
  126.       IK = (N*(N - 1))/2
  127.    60 IF (K .EQ. 0) GO TO 120
  128.          KK = IK + K
  129.          IKM1 = IK - (K - 1)
  130.          KS = 1
  131.          IF (KPVT(K) .LT. 0) KS = 2
  132.          KP = IABS(KPVT(K))
  133.          KPS = K + 1 - KS
  134.          IF (KP .EQ. KPS) GO TO 70
  135.             T = Z(KPS)
  136.             Z(KPS) = Z(KP)
  137.             Z(KP) = T
  138.    70    CONTINUE
  139.          IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,Z(K))
  140.          Z(K) = Z(K) + EK
  141.          CALL SAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1)
  142.          IF (KS .EQ. 1) GO TO 80
  143.             IF (Z(K-1) .NE. 0.0E0) EK = SIGN(EK,Z(K-1))
  144.             Z(K-1) = Z(K-1) + EK
  145.             CALL SAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1)
  146.    80    CONTINUE
  147.          IF (KS .EQ. 2) GO TO 100
  148.             IF (ABS(Z(K)) .LE. ABS(AP(KK))) GO TO 90
  149.                S = ABS(AP(KK))/ABS(Z(K))
  150.                CALL SSCAL(N,S,Z,1)
  151.                EK = S*EK
  152.    90       CONTINUE
  153.             IF (AP(KK) .NE. 0.0E0) Z(K) = Z(K)/AP(KK)
  154.             IF (AP(KK) .EQ. 0.0E0) Z(K) = 1.0E0
  155.          GO TO 110
  156.   100    CONTINUE
  157.             KM1K = IK + K - 1
  158.             KM1KM1 = IKM1 + K - 1
  159.             AK = AP(KK)/AP(KM1K)
  160.             AKM1 = AP(KM1KM1)/AP(KM1K)
  161.             BK = Z(K)/AP(KM1K)
  162.             BKM1 = Z(K-1)/AP(KM1K)
  163.             DENOM = AK*AKM1 - 1.0E0
  164.             Z(K) = (AKM1*BK - BKM1)/DENOM
  165.             Z(K-1) = (AK*BKM1 - BK)/DENOM
  166.   110    CONTINUE
  167.          K = K - KS
  168.          IK = IK - K
  169.          IF (KS .EQ. 2) IK = IK - (K + 1)
  170.       GO TO 60
  171.   120 CONTINUE
  172.       S = 1.0E0/SASUM(N,Z,1)
  173.       CALL SSCAL(N,S,Z,1)
  174. C
  175. C     SOLVE TRANS(U)*Y = W
  176. C
  177.       K = 1
  178.       IK = 0
  179.   130 IF (K .GT. N) GO TO 160
  180.          KS = 1
  181.          IF (KPVT(K) .LT. 0) KS = 2
  182.          IF (K .EQ. 1) GO TO 150
  183.             Z(K) = Z(K) + SDOT(K-1,AP(IK+1),1,Z(1),1)
  184.             IKP1 = IK + K
  185.             IF (KS .EQ. 2)
  186.      *         Z(K+1) = Z(K+1) + SDOT(K-1,AP(IKP1+1),1,Z(1),1)
  187.             KP = IABS(KPVT(K))
  188.             IF (KP .EQ. K) GO TO 140
  189.                T = Z(K)
  190.                Z(K) = Z(KP)
  191.                Z(KP) = T
  192.   140       CONTINUE
  193.   150    CONTINUE
  194.          IK = IK + K
  195.          IF (KS .EQ. 2) IK = IK + (K + 1)
  196.          K = K + KS
  197.       GO TO 130
  198.   160 CONTINUE
  199.       S = 1.0E0/SASUM(N,Z,1)
  200.       CALL SSCAL(N,S,Z,1)
  201. C
  202.       YNORM = 1.0E0
  203. C
  204. C     SOLVE U*D*V = Y
  205. C
  206.       K = N
  207.       IK = N*(N - 1)/2
  208.   170 IF (K .EQ. 0) GO TO 230
  209.          KK = IK + K
  210.          IKM1 = IK - (K - 1)
  211.          KS = 1
  212.          IF (KPVT(K) .LT. 0) KS = 2
  213.          IF (K .EQ. KS) GO TO 190
  214.             KP = IABS(KPVT(K))
  215.             KPS = K + 1 - KS
  216.             IF (KP .EQ. KPS) GO TO 180
  217.                T = Z(KPS)
  218.                Z(KPS) = Z(KP)
  219.                Z(KP) = T
  220.   180       CONTINUE
  221.             CALL SAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1)
  222.             IF (KS .EQ. 2) CALL SAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1)
  223.   190    CONTINUE
  224.          IF (KS .EQ. 2) GO TO 210
  225.             IF (ABS(Z(K)) .LE. ABS(AP(KK))) GO TO 200
  226.                S = ABS(AP(KK))/ABS(Z(K))
  227.                CALL SSCAL(N,S,Z,1)
  228.                YNORM = S*YNORM
  229.   200       CONTINUE
  230.             IF (AP(KK) .NE. 0.0E0) Z(K) = Z(K)/AP(KK)
  231.             IF (AP(KK) .EQ. 0.0E0) Z(K) = 1.0E0
  232.          GO TO 220
  233.   210    CONTINUE
  234.             KM1K = IK + K - 1
  235.             KM1KM1 = IKM1 + K - 1
  236.             AK = AP(KK)/AP(KM1K)
  237.             AKM1 = AP(KM1KM1)/AP(KM1K)
  238.             BK = Z(K)/AP(KM1K)
  239.             BKM1 = Z(K-1)/AP(KM1K)
  240.             DENOM = AK*AKM1 - 1.0E0
  241.             Z(K) = (AKM1*BK - BKM1)/DENOM
  242.             Z(K-1) = (AK*BKM1 - BK)/DENOM
  243.   220    CONTINUE
  244.          K = K - KS
  245.          IK = IK - K
  246.          IF (KS .EQ. 2) IK = IK - (K + 1)
  247.       GO TO 170
  248.   230 CONTINUE
  249.       S = 1.0E0/SASUM(N,Z,1)
  250.       CALL SSCAL(N,S,Z,1)
  251.       YNORM = S*YNORM
  252. C
  253. C     SOLVE TRANS(U)*Z = V
  254. C
  255.       K = 1
  256.       IK = 0
  257.   240 IF (K .GT. N) GO TO 270
  258.          KS = 1
  259.          IF (KPVT(K) .LT. 0) KS = 2
  260.          IF (K .EQ. 1) GO TO 260
  261.             Z(K) = Z(K) + SDOT(K-1,AP(IK+1),1,Z(1),1)
  262.             IKP1 = IK + K
  263.             IF (KS .EQ. 2)
  264.      *         Z(K+1) = Z(K+1) + SDOT(K-1,AP(IKP1+1),1,Z(1),1)
  265.             KP = IABS(KPVT(K))
  266.             IF (KP .EQ. K) GO TO 250
  267.                T = Z(K)
  268.                Z(K) = Z(KP)
  269.                Z(KP) = T
  270.   250       CONTINUE
  271.   260    CONTINUE
  272.          IK = IK + K
  273.          IF (KS .EQ. 2) IK = IK + (K + 1)
  274.          K = K + KS
  275.       GO TO 240
  276.   270 CONTINUE
  277. C     MAKE ZNORM = 1.0
  278.       S = 1.0E0/SASUM(N,Z,1)
  279.       CALL SSCAL(N,S,Z,1)
  280.       YNORM = S*YNORM
  281. C
  282.       IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
  283.       IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
  284.       RETURN
  285.       END
  286.