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 / WGEFA.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  2.8 KB  |  96 lines
  1.       SUBROUTINE WGEFA (AR, AI, LDA, N, IPVT, INFO)
  2.       IMPLICIT NONE
  3. C
  4. C FACTORS A DOUBLE-COMPLEX MATRIX BY GAUSSIAN ELIMINATION
  5. C  USUALLY CALLED BY WGECO, BUT IT CAN BE CALLED DIRECTLY WITH A SAVING IN
  6. C  TIME IF RCOND IS NOT NEEDED. (TIME FOR WGECO) = (1 + 9/N)*(TIME FOR WGEFA)
  7. C
  8. C ON ENTRY:
  9. C    A    DOUBLE-COMPLEX(LDA, N),  MATRIX TO BE FACTORED
  10. C    LDA  INTEGER,  LEADING DIMENSION OF THE ARRAY  A
  11. C    N    INTEGER,  ORDER OF THE MATRIX  A
  12. C
  13. C ON RETURN:
  14. C    A     AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS USED TO OBTAIN IT.
  15. C           THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE L  IS A PRODUCT
  16. C           OF PERMUTATION AND UNIT LOWER TRIANGULAR MATRICES AND  U  IS
  17. C           UPPER TRIANGULAR.
  18. C    IPVT  INTEGER(N),  AN INTEGER VECTOR OF PIVOT INDICES
  19. C    INFO  INTEGER
  20. C           = 0  NORMAL VALUE.
  21. C           = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR CONDITION FOR
  22. C                 THIS SUBROUTINE, BUT IT DOES INDICATE THAT WGESL OR WGEDI
  23. C                 WILL DIVIDE BY ZERO IF CALLED.  USE RCOND IN WGECO FOR A
  24. C                 RELIABLE INDICATION OF SINGULARITY
  25. C
  26. C VERSION 07/01/79 CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  27. C
  28.       INTEGER LDA, N, IPVT(*), INFO
  29.       DOUBLE PRECISION AR(LDA,*), AI(LDA,*)
  30. C
  31.       DOUBLE PRECISION TR, TI
  32.       INTEGER J, K, KP1, L, NM1
  33.       DOUBLE PRECISION ZDUMR, ZDUMI
  34. C
  35.       INTEGER IWAMAX
  36.       DOUBLE PRECISION CABS1
  37. C
  38. C
  39.       CABS1 (ZDUMR, ZDUMI) = DABS (ZDUMR)+DABS (ZDUMI)
  40. C
  41. C
  42. C ***      GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
  43.       INFO = 0
  44.       NM1 = N-1
  45.       IF (NM1.LT.1) GO TO 70
  46.       DO 60 K = 1, NM1
  47.         KP1 = K+1
  48. C
  49. C ***      FIND L = PIVOT INDEX
  50.         L = IWAMAX (N-K+1, AR(K,K), AI(K,K), 1)+K-1
  51.         IPVT(K) = L
  52. C
  53. C ***      ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
  54.         IF (CABS1 (AR(L,K), AI(L,K)).EQ.0.0D0) GO TO 40
  55. C
  56. C ***      INTERCHANGE IF NECESSARY
  57.         IF (L.EQ.K) GO TO 10
  58.         TR = AR(L,K)
  59.         TI = AI(L,K)
  60.         AR(L,K) = AR(K,K)
  61.         AI(L,K) = AI(K,K)
  62.         AR(K,K) = TR
  63.         AI(K,K) = TI
  64. 10      CONTINUE
  65. C
  66. C ***      COMPUTE MULTIPLIERS
  67.         CALL WDIV (-1.0D0, 0.0D0, AR(K,K), AI(K,K), TR, TI)
  68.         CALL WSCAL (N-K, TR, TI, AR(K+1,K), AI(K+1,K), 1)
  69. C
  70. C ***      ROW ELIMINATION WITH COLUMN INDEXING
  71.         DO 30 J = KP1, N
  72.           TR = AR(L,J)
  73.           TI = AI(L,J)
  74.           IF (L.EQ.K) GO TO 20
  75.           AR(L,J) = AR(K,J)
  76.           AI(L,J) = AI(K,J)
  77.           AR(K,J) = TR
  78.           AI(K,J) = TI
  79. 20        CONTINUE
  80.           CALL WAXPY (N-K, TR, TI, AR(K+1,K), AI(K+1,K),
  81.      .                1, AR(K+1,J), AI(K+1,J), 1)
  82. 30      CONTINUE
  83.         GO TO 50
  84. C
  85. 40      CONTINUE
  86.         INFO = K
  87. 50      CONTINUE
  88. 60    CONTINUE
  89. C
  90. 70    CONTINUE
  91.       IPVT(N) = N
  92.       IF (CABS1 (AR(N,N), AI(N,N)).EQ.0.0D0) INFO = N
  93. C
  94.       RETURN
  95.       END
  96.