home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / linpack / dgbfa.f next >
Text File  |  1996-09-28  |  5KB  |  175 lines
  1.       SUBROUTINE DGBFA(ABD,LDA,N,ML,MU,IPVT,INFO)
  2.       INTEGER LDA,N,ML,MU,IPVT(1),INFO
  3.       DOUBLE PRECISION ABD(LDA,1)
  4. C
  5. C     DGBFA FACTORS A DOUBLE PRECISION BAND MATRIX BY ELIMINATION.
  6. C
  7. C     DGBFA IS USUALLY CALLED BY DGBCO, BUT IT CAN BE CALLED
  8. C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
  9. C
  10. C     ON ENTRY
  11. C
  12. C        ABD     DOUBLE PRECISION(LDA, N)
  13. C                CONTAINS THE MATRIX IN BAND STORAGE.  THE COLUMNS
  14. C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  ABD  AND
  15. C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
  16. C                ML+1 THROUGH 2*ML+MU+1 OF  ABD .
  17. C                SEE THE COMMENTS BELOW FOR DETAILS.
  18. C
  19. C        LDA     INTEGER
  20. C                THE LEADING DIMENSION OF THE ARRAY  ABD .
  21. C                LDA MUST BE .GE. 2*ML + MU + 1 .
  22. C
  23. C        N       INTEGER
  24. C                THE ORDER OF THE ORIGINAL MATRIX.
  25. C
  26. C        ML      INTEGER
  27. C                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
  28. C                0 .LE. ML .LT. N .
  29. C
  30. C        MU      INTEGER
  31. C                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
  32. C                0 .LE. MU .LT. N .
  33. C                MORE EFFICIENT IF  ML .LE. MU .
  34. C     ON RETURN
  35. C
  36. C        ABD     AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
  37. C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
  38. C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
  39. C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
  40. C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
  41. C
  42. C        IPVT    INTEGER(N)
  43. C                AN INTEGER VECTOR OF PIVOT INDICES.
  44. C
  45. C        INFO    INTEGER
  46. C                = 0  NORMAL VALUE.
  47. C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
  48. C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
  49. C                     INDICATE THAT DGBSL WILL DIVIDE BY ZERO IF
  50. C                     CALLED.  USE  RCOND  IN DGBCO FOR A RELIABLE
  51. C                     INDICATION OF SINGULARITY.
  52. C
  53. C     BAND STORAGE
  54. C
  55. C           IF  A  IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT
  56. C           WILL SET UP THE INPUT.
  57. C
  58. C                   ML = (BAND WIDTH BELOW THE DIAGONAL)
  59. C                   MU = (BAND WIDTH ABOVE THE DIAGONAL)
  60. C                   M = ML + MU + 1
  61. C                   DO 20 J = 1, N
  62. C                      I1 = MAX0(1, J-MU)
  63. C                      I2 = MIN0(N, J+ML)
  64. C                      DO 10 I = I1, I2
  65. C                         K = I - J + M
  66. C                         ABD(K,J) = A(I,J)
  67. C                10    CONTINUE
  68. C                20 CONTINUE
  69. C
  70. C           THIS USES ROWS  ML+1  THROUGH  2*ML+MU+1  OF  ABD .
  71. C           IN ADDITION, THE FIRST  ML  ROWS IN  ABD  ARE USED FOR
  72. C           ELEMENTS GENERATED DURING THE TRIANGULARIZATION.
  73. C           THE TOTAL NUMBER OF ROWS NEEDED IN  ABD  IS  2*ML+MU+1 .
  74. C           THE  ML+MU BY ML+MU  UPPER LEFT TRIANGLE AND THE
  75. C           ML BY ML  LOWER RIGHT TRIANGLE ARE NOT REFERENCED.
  76. C
  77. C     LINPACK. THIS VERSION DATED 08/14/78 .
  78. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  79. C
  80. C     SUBROUTINES AND FUNCTIONS
  81. C
  82. C     BLAS DAXPY,DSCAL,IDAMAX
  83. C     FORTRAN MAX0,MIN0
  84. C
  85. C     INTERNAL VARIABLES
  86. C
  87.       DOUBLE PRECISION T
  88.       INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
  89. C
  90. C
  91.       M = ML + MU + 1
  92.       INFO = 0
  93. C
  94. C     ZERO INITIAL FILL-IN COLUMNS
  95. C
  96.       J0 = MU + 2
  97.       J1 = MIN0(N,M) - 1
  98.       IF (J1 .LT. J0) GO TO 30
  99.       DO 20 JZ = J0, J1
  100.          I0 = M + 1 - JZ
  101.          DO 10 I = I0, ML
  102.             ABD(I,JZ) = 0.0D0
  103.    10    CONTINUE
  104.    20 CONTINUE
  105.    30 CONTINUE
  106.       JZ = J1
  107.       JU = 0
  108. C
  109. C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
  110. C
  111.       NM1 = N - 1
  112.       IF (NM1 .LT. 1) GO TO 130
  113.       DO 120 K = 1, NM1
  114.          KP1 = K + 1
  115. C
  116. C        ZERO NEXT FILL-IN COLUMN
  117. C
  118.          JZ = JZ + 1
  119.          IF (JZ .GT. N) GO TO 50
  120.          IF (ML .LT. 1) GO TO 50
  121.             DO 40 I = 1, ML
  122.                ABD(I,JZ) = 0.0D0
  123.    40       CONTINUE
  124.    50    CONTINUE
  125. C
  126. C        FIND L = PIVOT INDEX
  127. C
  128.          LM = MIN0(ML,N-K)
  129.          L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
  130.          IPVT(K) = L + K - M
  131. C
  132. C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
  133. C
  134.          IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
  135. C
  136. C           INTERCHANGE IF NECESSARY
  137. C
  138.             IF (L .EQ. M) GO TO 60
  139.                T = ABD(L,K)
  140.                ABD(L,K) = ABD(M,K)
  141.                ABD(M,K) = T
  142.    60       CONTINUE
  143. C
  144. C           COMPUTE MULTIPLIERS
  145. C
  146.             T = -1.0D0/ABD(M,K)
  147.             CALL DSCAL(LM,T,ABD(M+1,K),1)
  148. C
  149. C           ROW ELIMINATION WITH COLUMN INDEXING
  150. C
  151.             JU = MIN0(MAX0(JU,MU+IPVT(K)),N)
  152.             MM = M
  153.             IF (JU .LT. KP1) GO TO 90
  154.             DO 80 J = KP1, JU
  155.                L = L - 1
  156.                MM = MM - 1
  157.                T = ABD(L,J)
  158.                IF (L .EQ. MM) GO TO 70
  159.                   ABD(L,J) = ABD(MM,J)
  160.                   ABD(MM,J) = T
  161.    70          CONTINUE
  162.                CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
  163.    80       CONTINUE
  164.    90       CONTINUE
  165.          GO TO 110
  166.   100    CONTINUE
  167.             INFO = K
  168.   110    CONTINUE
  169.   120 CONTINUE
  170.   130 CONTINUE
  171.       IPVT(N) = N
  172.       IF (ABD(M,N) .EQ. 0.0D0) INFO = N
  173.       RETURN
  174.       END
  175.