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 / dgbsl.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  136 lines
  1.       SUBROUTINE DGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB)
  2.       INTEGER LDA,N,ML,MU,IPVT(1),JOB
  3.       DOUBLE PRECISION ABD(LDA,1),B(1)
  4. C
  5. C     DGBSL SOLVES THE DOUBLE PRECISION BAND SYSTEM
  6. C     A * X = B  OR  TRANS(A) * X = B
  7. C     USING THE FACTORS COMPUTED BY DGBCO OR DGBFA.
  8. C
  9. C     ON ENTRY
  10. C
  11. C        ABD     DOUBLE PRECISION(LDA, N)
  12. C                THE OUTPUT FROM DGBCO OR DGBFA.
  13. C
  14. C        LDA     INTEGER
  15. C                THE LEADING DIMENSION OF THE ARRAY  ABD .
  16. C
  17. C        N       INTEGER
  18. C                THE ORDER OF THE ORIGINAL MATRIX.
  19. C
  20. C        ML      INTEGER
  21. C                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
  22. C
  23. C        MU      INTEGER
  24. C                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
  25. C
  26. C        IPVT    INTEGER(N)
  27. C                THE PIVOT VECTOR FROM DGBCO OR DGBFA.
  28. C
  29. C        B       DOUBLE PRECISION(N)
  30. C                THE RIGHT HAND SIDE VECTOR.
  31. C
  32. C        JOB     INTEGER
  33. C                = 0         TO SOLVE  A*X = B ,
  34. C                = NONZERO   TO SOLVE  TRANS(A)*X = B , WHERE
  35. C                            TRANS(A)  IS THE TRANSPOSE.
  36. C
  37. C     ON RETURN
  38. C
  39. C        B       THE SOLUTION VECTOR  X .
  40. C
  41. C     ERROR CONDITION
  42. C
  43. C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
  44. C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
  45. C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
  46. C        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE
  47. C        CALLED CORRECTLY AND IF DGBCO HAS SET RCOND .GT. 0.0
  48. C        OR DGBFA HAS SET INFO .EQ. 0 .
  49. C
  50. C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
  51. C     WITH  P  COLUMNS
  52. C           CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
  53. C           IF (RCOND IS TOO SMALL) GO TO ...
  54. C           DO 10 J = 1, P
  55. C              CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
  56. C        10 CONTINUE
  57. C
  58. C     LINPACK. THIS VERSION DATED 08/14/78 .
  59. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  60. C
  61. C     SUBROUTINES AND FUNCTIONS
  62. C
  63. C     BLAS DAXPY,DDOT
  64. C     FORTRAN MIN0
  65. C
  66. C     INTERNAL VARIABLES
  67. C
  68.       DOUBLE PRECISION DDOT,T
  69.       INTEGER K,KB,L,LA,LB,LM,M,NM1
  70. C
  71.       M = MU + ML + 1
  72.       NM1 = N - 1
  73.       IF (JOB .NE. 0) GO TO 50
  74. C
  75. C        JOB = 0 , SOLVE  A * X = B
  76. C        FIRST SOLVE L*Y = B
  77. C
  78.          IF (ML .EQ. 0) GO TO 30
  79.          IF (NM1 .LT. 1) GO TO 30
  80.             DO 20 K = 1, NM1
  81.                LM = MIN0(ML,N-K)
  82.                L = IPVT(K)
  83.                T = B(L)
  84.                IF (L .EQ. K) GO TO 10
  85.                   B(L) = B(K)
  86.                   B(K) = T
  87.    10          CONTINUE
  88.                CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
  89.    20       CONTINUE
  90.    30    CONTINUE
  91. C
  92. C        NOW SOLVE  U*X = Y
  93. C
  94.          DO 40 KB = 1, N
  95.             K = N + 1 - KB
  96.             B(K) = B(K)/ABD(M,K)
  97.             LM = MIN0(K,M) - 1
  98.             LA = M - LM
  99.             LB = K - LM
  100.             T = -B(K)
  101.             CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
  102.    40    CONTINUE
  103.       GO TO 100
  104.    50 CONTINUE
  105. C
  106. C        JOB = NONZERO, SOLVE  TRANS(A) * X = B
  107. C        FIRST SOLVE  TRANS(U)*Y = B
  108. C
  109.          DO 60 K = 1, N
  110.             LM = MIN0(K,M) - 1
  111.             LA = M - LM
  112.             LB = K - LM
  113.             T = DDOT(LM,ABD(LA,K),1,B(LB),1)
  114.             B(K) = (B(K) - T)/ABD(M,K)
  115.    60    CONTINUE
  116. C
  117. C        NOW SOLVE TRANS(L)*X = Y
  118. C
  119.          IF (ML .EQ. 0) GO TO 90
  120.          IF (NM1 .LT. 1) GO TO 90
  121.             DO 80 KB = 1, NM1
  122.                K = N - KB
  123.                LM = MIN0(ML,N-K)
  124.                B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
  125.                L = IPVT(K)
  126.                IF (L .EQ. K) GO TO 70
  127.                   T = B(L)
  128.                   B(L) = B(K)
  129.                   B(K) = T
  130.    70          CONTINUE
  131.    80       CONTINUE
  132.    90    CONTINUE
  133.   100 CONTINUE
  134.       RETURN
  135.       END
  136.