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