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