home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / fortran / linpklib.arc / SPPCO.FOR < prev    next >
Text File  |  1984-01-06  |  7KB  |  220 lines
  1.       SUBROUTINE SPPCO(AP,N,RCOND,Z,INFO)
  2.       INTEGER N,INFO
  3.       REAL AP(1),Z(1)
  4.       REAL RCOND
  5. C
  6. C     SPPCO FACTORS A REAL SYMMETRIC POSITIVE DEFINITE
  7. C     MATRIX STORED IN PACKED FORM
  8. C     AND ESTIMATES THE CONDITION OF THE MATRIX.
  9. C
  10. C     IF  RCOND  IS NOT NEEDED, SPPFA IS SLIGHTLY FASTER.
  11. C     TO SOLVE  A*X = B , FOLLOW SPPCO BY SPPSL.
  12. C     TO COMPUTE  INVERSE(A)*C , FOLLOW SPPCO BY SPPSL.
  13. C     TO COMPUTE  DETERMINANT(A) , FOLLOW SPPCO BY SPPDI.
  14. C     TO COMPUTE  INVERSE(A) , FOLLOW SPPCO BY SPPDI.
  15. C
  16. C     ON ENTRY
  17. C
  18. C        AP      REAL (N*(N+1)/2)
  19. C                THE PACKED FORM OF A SYMMETRIC MATRIX  A .  THE
  20. C                COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY
  21. C                IN A ONE-DIMENSIONAL ARRAY OF LENGTH  N*(N+1)/2 .
  22. C                SEE COMMENTS BELOW FOR DETAILS.
  23. C
  24. C        N       INTEGER
  25. C                THE ORDER OF THE MATRIX  A .
  26. C
  27. C     ON RETURN
  28. C
  29. C        AP      AN UPPER TRIANGULAR MATRIX  R , STORED IN PACKED
  30. C                FORM, SO THAT  A = TRANS(R)*R .
  31. C                IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
  32. C
  33. C        RCOND   REAL
  34. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .
  35. C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS
  36. C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE
  37. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
  38. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
  39. C                           1.0 + RCOND .EQ. 1.0
  40. C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING
  41. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
  42. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
  43. C                UNDERFLOWS.  IF INFO .NE. 0 , RCOND IS UNCHANGED.
  44. C
  45. C        Z       REAL(N)
  46. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
  47. C                IF  A  IS SINGULAR TO WORKING PRECISION, THEN  Z  IS
  48. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
  49. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
  50. C                IF  INFO .NE. 0 , Z  IS UNCHANGED.
  51. C
  52. C        INFO    INTEGER
  53. C                = 0  FOR NORMAL RETURN.
  54. C                = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
  55. C                     OF ORDER  K  IS NOT POSITIVE DEFINITE.
  56. C
  57. C     PACKED STORAGE
  58. C
  59. C          THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER
  60. C          TRIANGLE OF A SYMMETRIC MATRIX.
  61. C
  62. C                K = 0
  63. C                DO 20 J = 1, N
  64. C                   DO 10 I = 1, J
  65. C                      K = K + 1
  66. C                      AP(K) = A(I,J)
  67. C             10    CONTINUE
  68. C             20 CONTINUE
  69. C
  70. C     LINPACK.  THIS VERSION DATED 08/14/78 .
  71. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  72. C
  73. C     SUBROUTINES AND FUNCTIONS
  74. C
  75. C     LINPACK SPPFA
  76. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  77. C     FORTRAN ABS,AMAX1,REAL,SIGN
  78. C
  79. C     INTERNAL VARIABLES
  80. C
  81.       REAL SDOT,EK,T,WK,WKM
  82.       REAL ANORM,S,SASUM,SM,YNORM
  83.       INTEGER I,IJ,J,JM1,J1,K,KB,KJ,KK,KP1
  84. C
  85. C
  86. C     FIND NORM OF A
  87. C
  88.       J1 = 1
  89.       DO 30 J = 1, N
  90.          Z(J) = SASUM(J,AP(J1),1)
  91.          IJ = J1
  92.          J1 = J1 + J
  93.          JM1 = J - 1
  94.          IF (JM1 .LT. 1) GO TO 20
  95.          DO 10 I = 1, JM1
  96.             Z(I) = Z(I) + ABS(AP(IJ))
  97.             IJ = IJ + 1
  98.    10    CONTINUE
  99.    20    CONTINUE
  100.    30 CONTINUE
  101.       ANORM = 0.0E0
  102.       DO 40 J = 1, N
  103.          ANORM = AMAX1(ANORM,Z(J))
  104.    40 CONTINUE
  105. C
  106. C     FACTOR
  107. C
  108.       CALL SPPFA(AP,N,INFO)
  109.       IF (INFO .NE. 0) GO TO 180
  110. C
  111. C        RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
  112. C        ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
  113. C        THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
  114. C        GROWTH IN THE ELEMENTS OF W  WHERE  TRANS(R)*W = E .
  115. C        THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
  116. C
  117. C        SOLVE TRANS(R)*W = E
  118. C
  119.          EK = 1.0E0
  120.          DO 50 J = 1, N
  121.             Z(J) = 0.0E0
  122.    50    CONTINUE
  123.          KK = 0
  124.          DO 110 K = 1, N
  125.             KK = KK + K
  126.             IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K))
  127.             IF (ABS(EK-Z(K)) .LE. AP(KK)) GO TO 60
  128.                S = AP(KK)/ABS(EK-Z(K))
  129.                CALL SSCAL(N,S,Z,1)
  130.                EK = S*EK
  131.    60       CONTINUE
  132.             WK = EK - Z(K)
  133.             WKM = -EK - Z(K)
  134.             S = ABS(WK)
  135.             SM = ABS(WKM)
  136.             WK = WK/AP(KK)
  137.             WKM = WKM/AP(KK)
  138.             KP1 = K + 1
  139.             KJ = KK + K
  140.             IF (KP1 .GT. N) GO TO 100
  141.                DO 70 J = KP1, N
  142.                   SM = SM + ABS(Z(J)+WKM*AP(KJ))
  143.                   Z(J) = Z(J) + WK*AP(KJ)
  144.                   S = S + ABS(Z(J))
  145.                   KJ = KJ + J
  146.    70          CONTINUE
  147.                IF (S .GE. SM) GO TO 90
  148.                   T = WKM - WK
  149.                   WK = WKM
  150.                   KJ = KK + K
  151.                   DO 80 J = KP1, N
  152.                      Z(J) = Z(J) + T*AP(KJ)
  153.                      KJ = KJ + J
  154.    80             CONTINUE
  155.    90          CONTINUE
  156.   100       CONTINUE
  157.             Z(K) = WK
  158.   110    CONTINUE
  159.          S = 1.0E0/SASUM(N,Z,1)
  160.          CALL SSCAL(N,S,Z,1)
  161. C
  162. C        SOLVE R*Y = W
  163. C
  164.          DO 130 KB = 1, N
  165.             K = N + 1 - KB
  166.             IF (ABS(Z(K)) .LE. AP(KK)) GO TO 120
  167.                S = AP(KK)/ABS(Z(K))
  168.                CALL SSCAL(N,S,Z,1)
  169.   120       CONTINUE
  170.             Z(K) = Z(K)/AP(KK)
  171.             KK = KK - K
  172.             T = -Z(K)
  173.             CALL SAXPY(K-1,T,AP(KK+1),1,Z(1),1)
  174.   130    CONTINUE
  175.          S = 1.0E0/SASUM(N,Z,1)
  176.          CALL SSCAL(N,S,Z,1)
  177. C
  178.          YNORM = 1.0E0
  179. C
  180. C        SOLVE TRANS(R)*V = Y
  181. C
  182.          DO 150 K = 1, N
  183.             Z(K) = Z(K) - SDOT(K-1,AP(KK+1),1,Z(1),1)
  184.             KK = KK + K
  185.             IF (ABS(Z(K)) .LE. AP(KK)) GO TO 140
  186.                S = AP(KK)/ABS(Z(K))
  187.                CALL SSCAL(N,S,Z,1)
  188.                YNORM = S*YNORM
  189.   140       CONTINUE
  190.             Z(K) = Z(K)/AP(KK)
  191.   150    CONTINUE
  192.          S = 1.0E0/SASUM(N,Z,1)
  193.          CALL SSCAL(N,S,Z,1)
  194.          YNORM = S*YNORM
  195. C
  196. C        SOLVE R*Z = V
  197. C
  198.          DO 170 KB = 1, N
  199.             K = N + 1 - KB
  200.             IF (ABS(Z(K)) .LE. AP(KK)) GO TO 160
  201.                S = AP(KK)/ABS(Z(K))
  202.                CALL SSCAL(N,S,Z,1)
  203.                YNORM = S*YNORM
  204.   160       CONTINUE
  205.             Z(K) = Z(K)/AP(KK)
  206.             KK = KK - K
  207.             T = -Z(K)
  208.             CALL SAXPY(K-1,T,AP(KK+1),1,Z(1),1)
  209.   170    CONTINUE
  210. C        MAKE ZNORM = 1.0
  211.          S = 1.0E0/SASUM(N,Z,1)
  212.          CALL SSCAL(N,S,Z,1)
  213.          YNORM = S*YNORM
  214. C
  215.          IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
  216.          IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
  217.   180 CONTINUE
  218.       RETURN
  219.       END
  220.