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 / dgedi.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  129 lines
  1.       SUBROUTINE DGEDI(A,LDA,N,IPVT,DET,WORK,JOB)
  2.       INTEGER LDA,N,IPVT(1),JOB
  3.       DOUBLE PRECISION A(LDA,1),DET(2),WORK(1)
  4. C
  5. C     DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX
  6. C     USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
  7. C
  8. C     ON ENTRY
  9. C
  10. C        A       DOUBLE PRECISION(LDA, N)
  11. C                THE OUTPUT FROM DGECO OR DGEFA.
  12. C
  13. C        LDA     INTEGER
  14. C                THE LEADING DIMENSION OF THE ARRAY  A .
  15. C
  16. C        N       INTEGER
  17. C                THE ORDER OF THE MATRIX  A .
  18. C
  19. C        IPVT    INTEGER(N)
  20. C                THE PIVOT VECTOR FROM DGECO OR DGEFA.
  21. C
  22. C        WORK    DOUBLE PRECISION(N)
  23. C                WORK VECTOR.  CONTENTS DESTROYED.
  24. C
  25. C        JOB     INTEGER
  26. C                = 11   BOTH DETERMINANT AND INVERSE.
  27. C                = 01   INVERSE ONLY.
  28. C                = 10   DETERMINANT ONLY.
  29. C
  30. C     ON RETURN
  31. C
  32. C        A       INVERSE OF ORIGINAL MATRIX IF REQUESTED.
  33. C                OTHERWISE UNCHANGED.
  34. C
  35. C        DET     DOUBLE PRECISION(2)
  36. C                DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.
  37. C                OTHERWISE NOT REFERENCED.
  38. C                DETERMINANT = DET(1) * 10.0**DET(2)
  39. C                WITH  1.0 .LE. DABS(DET(1)) .LT. 10.0
  40. C                OR  DET(1) .EQ. 0.0 .
  41. C
  42. C     ERROR CONDITION
  43. C
  44. C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
  45. C        A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.
  46. C        IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY
  47. C        AND IF DGECO HAS SET RCOND .GT. 0.0 OR DGEFA HAS SET
  48. C        INFO .EQ. 0 .
  49. C
  50. C     LINPACK. THIS VERSION DATED 08/14/78 .
  51. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  52. C
  53. C     SUBROUTINES AND FUNCTIONS
  54. C
  55. C     BLAS DAXPY,DSCAL,DSWAP
  56. C     FORTRAN DABS,MOD
  57. C
  58. C     INTERNAL VARIABLES
  59. C
  60.       DOUBLE PRECISION T
  61.       DOUBLE PRECISION TEN
  62.       INTEGER I,J,K,KB,KP1,L,NM1
  63. C
  64. C
  65. C     COMPUTE DETERMINANT
  66. C
  67.       IF (JOB/10 .EQ. 0) GO TO 70
  68.          DET(1) = 1.0D0
  69.          DET(2) = 0.0D0
  70.          TEN = 10.0D0
  71.          DO 50 I = 1, N
  72.             IF (IPVT(I) .NE. I) DET(1) = -DET(1)
  73.             DET(1) = A(I,I)*DET(1)
  74. C        ...EXIT
  75.             IF (DET(1) .EQ. 0.0D0) GO TO 60
  76.    10       IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20
  77.                DET(1) = TEN*DET(1)
  78.                DET(2) = DET(2) - 1.0D0
  79.             GO TO 10
  80.    20       CONTINUE
  81.    30       IF (DABS(DET(1)) .LT. TEN) GO TO 40
  82.                DET(1) = DET(1)/TEN
  83.                DET(2) = DET(2) + 1.0D0
  84.             GO TO 30
  85.    40       CONTINUE
  86.    50    CONTINUE
  87.    60    CONTINUE
  88.    70 CONTINUE
  89. C
  90. C     COMPUTE INVERSE(U)
  91. C
  92.       IF (MOD(JOB,10) .EQ. 0) GO TO 150
  93.          DO 100 K = 1, N
  94.             A(K,K) = 1.0D0/A(K,K)
  95.             T = -A(K,K)
  96.             CALL DSCAL(K-1,T,A(1,K),1)
  97.             KP1 = K + 1
  98.             IF (N .LT. KP1) GO TO 90
  99.             DO 80 J = KP1, N
  100.                T = A(K,J)
  101.                A(K,J) = 0.0D0
  102.                CALL DAXPY(K,T,A(1,K),1,A(1,J),1)
  103.    80       CONTINUE
  104.    90       CONTINUE
  105.   100    CONTINUE
  106. C
  107. C        FORM INVERSE(U)*INVERSE(L)
  108. C
  109.          NM1 = N - 1
  110.          IF (NM1 .LT. 1) GO TO 140
  111.          DO 130 KB = 1, NM1
  112.             K = N - KB
  113.             KP1 = K + 1
  114.             DO 110 I = KP1, N
  115.                WORK(I) = A(I,K)
  116.                A(I,K) = 0.0D0
  117.   110       CONTINUE
  118.             DO 120 J = KP1, N
  119.                T = WORK(J)
  120.                CALL DAXPY(N,T,A(1,J),1,A(1,K),1)
  121.   120       CONTINUE
  122.             L = IPVT(K)
  123.             IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1)
  124.   130    CONTINUE
  125.   140    CONTINUE
  126.   150 CONTINUE
  127.       RETURN
  128.       END
  129.