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