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