home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / libs / matlab.lzh / MATLAB / MATLAB.LZH / Source / MatLab / WGEDI.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  3.8 KB  |  133 lines
  1.       SUBROUTINE WGEDI (AR, AI, LDA, N, IPVT, DETR, DETI,
  2.      .                  WORKR, WORKI, JOB)
  3.       IMPLICIT NONE
  4. C
  5. C COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX USING THE FACTORS
  6. C  COMPUTED BY WGECO OR WGEFA.
  7. C
  8. C ON ENTRY:
  9. C    A     DOUBLE-COMPLEX(LDA, N),  OUTPUT FROM WGECO OR WGEFA
  10. C    LDA   INTEGER,  LEADING DIMENSION OF THE ARRAY A
  11. C    N     INTEGER,  ORDER OF THE MATRIX A
  12. C    IPVT  INTEGER(N),  PIVOT VECTOR FROM WGECO OR WGEFA
  13. C    WORK  DOUBLE-COMPLEX(N),  WORK VECTOR.  CONTENTS DESTROYED
  14. C    JOB   INTEGER
  15. C            = 11  BOTH DETERMINANT AND INVERSE
  16. C            = 01  INVERSE ONLY
  17. C            = 10  DETERMINANT ONLY
  18. C
  19. C ON RETURN:
  20. C    A     INVERSE OF ORIGINAL MATRIX IF REQUESTED.  OTHERWISE UNCHANGED
  21. C
  22. C    DET   DOUBLE-COMPLEX(2),  DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.
  23. C           OTHERWISE NOT REFERENCED.  DETERMINANT = DET(1) * 10.0**DET(2)
  24. C           WITH  1.0 .LE. CABS1(DET(1) .LT. 10.0  OR  DET(1) .EQ. 0.0
  25. C
  26. C ERROR CONDITION
  27. C    DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A ZERO ON THE
  28. C    DIAGONAL AND THE INVERSE IS REQUESTED.  IT WILL NOT OCCUR IF THE
  29. C    SUBROUTINES ARE CALLED CORRECTLY AND IF WGECO HAS SET RCOND .GT. 0.0
  30. C    OR WGEFA HAS SET INFO .EQ. 0
  31. C
  32. C VERSION 07/01/79 CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  33. C
  34.       INTEGER LDA, N, IPVT(*), JOB
  35.       DOUBLE PRECISION AR(LDA,*), AI(LDA,*), DETR(*), DETI(*),
  36.      .                 WORKR(*), WORKI(*)
  37. C
  38.       DOUBLE PRECISION TR, TI, TEN
  39.       INTEGER I, J, K, KB, KP1, L, NM1
  40.       DOUBLE PRECISION ZDUMR, ZDUMI
  41. C
  42.       DOUBLE PRECISION CABS1
  43. C
  44. C
  45.       CABS1 (ZDUMR, ZDUMI) = DABS (ZDUMR)+DABS (ZDUMI)
  46. C
  47. C
  48. C ***      COMPUTE DETERMINANT
  49.       IF (JOB/10.EQ.0) GO TO 80
  50.       DETR(1) = 1.0D0
  51.       DETI(1) = 0.0D0
  52.       DETR(2) = 0.0D0
  53.       DETI(2) = 0.0D0
  54.       TEN = 10.0D0
  55.       DO 60 I = 1, N
  56.         IF (IPVT(I).EQ.I) GO TO 10
  57.         DETR(1) = -DETR(1)
  58.         DETI(1) = -DETI(1)
  59. 10      CONTINUE
  60.         CALL WMUL (AR(I,I), AI(I,I), DETR(1), DETI(1),
  61.      .             DETR(1), DETI(1))
  62. C
  63.         IF (CABS1 (DETR(1), DETI(1)).EQ.0.0D0) GO TO 70
  64. 20      CONTINUE
  65.         IF (CABS1 (DETR(1), DETI(1)).GE.1.0D0) GO TO 30
  66.         DETR(1) = TEN*DETR(1)
  67.         DETI(1) = TEN*DETI(1)
  68.         DETR(2) = DETR(2)-1.0D0
  69.         DETI(2) = DETI(2)-0.0D0
  70.         GO TO 20
  71. C
  72. 30      CONTINUE
  73. 40      CONTINUE
  74.         IF (CABS1 (DETR(1), DETI(1)).LT.TEN) GO TO 50
  75.         DETR(1) = DETR(1)/TEN
  76.         DETI(1) = DETI(1)/TEN
  77.         DETR(2) = DETR(2)+1.0D0
  78.         DETI(2) = DETI(2)+0.0D0
  79.         GO TO 40
  80. C
  81. 50      CONTINUE
  82. 60    CONTINUE
  83. 70    CONTINUE
  84. 80    CONTINUE
  85. C
  86. C ***      COMPUTE INVERSE(U)
  87.       IF (MOD (JOB, 10).EQ.0) GO TO 160
  88.       DO 110 K = 1, N
  89.         CALL WDIV (1.0D0, 0.0D0, AR(K,K), AI(K,K), AR(K,K), AI(K,K))
  90.         TR = -AR(K,K)
  91.         TI = -AI(K,K)
  92.         CALL WSCAL (K-1, TR, TI, AR(1,K), AI(1,K), 1)
  93.         KP1 = K+1
  94.         IF (N.LT.KP1) GO TO 100
  95.         DO 90 J = KP1, N
  96.           TR = AR(K,J)
  97.           TI = AI(K,J)
  98.           AR(K,J) = 0.0D0
  99.           AI(K,J) = 0.0D0
  100.           CALL WAXPY (K, TR, TI, AR(1,K), AI(1,K),
  101.      .                1, AR(1,J), AI(1,J), 1)
  102. 90      CONTINUE
  103. 100     CONTINUE
  104. 110   CONTINUE
  105. C
  106. C ***      FORM INVERSE(U)*INVERSE(L)
  107.       NM1 = N-1
  108.       IF (NM1.LT.1) GO TO 150
  109.       DO 140 KB = 1, NM1
  110.         K = N-KB
  111.         KP1 = K+1
  112.         DO 120 I = KP1, N
  113.           WORKR(I) = AR(I,K)
  114.           WORKI(I) = AI(I,K)
  115.           AR(I,K) = 0.0D0
  116.           AI(I,K) = 0.0D0
  117. 120     CONTINUE
  118.         DO 130 J = KP1, N
  119.           TR = WORKR(J)
  120.           TI = WORKI(J)
  121.           CALL WAXPY (N, TR, TI, AR(1,J), AI(1,J),
  122.      .                1, AR(1,K), AI(1,K), 1)
  123. 130     CONTINUE
  124.         L = IPVT(K)
  125.         IF (L.NE.K) CALL WSWAP (N, AR(1,K), AI(1,K), 1,
  126.      .                             AR(1,L), AI(1,L), 1)
  127. 140   CONTINUE
  128. 150   CONTINUE
  129. 160   CONTINUE
  130. C
  131.       RETURN
  132.       END
  133.