cgefa.F 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108
  1. SUBROUTINE CGEFA(A,LDA,N,IPVT,INFO)
  2. INTEGER LDA,N,IPVT(*),INFO
  3. COMPLEX A(LDA,*)
  4. C
  5. C CGEFA FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION.
  6. C
  7. C CGEFA IS USUALLY CALLED BY CGECO, BUT IT CAN BE CALLED
  8. C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
  9. C (TIME FOR CGECO) = (1 + 9/N)*(TIME FOR CGEFA) .
  10. C
  11. C ON ENTRY
  12. C
  13. C A COMPLEX(LDA, N)
  14. C THE MATRIX TO BE FACTORED.
  15. C
  16. C LDA INTEGER
  17. C THE LEADING DIMENSION OF THE ARRAY A .
  18. C
  19. C N INTEGER
  20. C THE ORDER OF THE MATRIX A .
  21. C
  22. C ON RETURN
  23. C
  24. C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
  25. C WHICH WERE USED TO OBTAIN IT.
  26. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
  27. C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
  28. C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
  29. C
  30. C IPVT INTEGER(N)
  31. C AN INTEGER VECTOR OF PIVOT INDICES.
  32. C
  33. C INFO INTEGER
  34. C = 0 NORMAL VALUE.
  35. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
  36. C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
  37. C INDICATE THAT CGESL OR CGEDI WILL DIVIDE BY ZERO
  38. C IF CALLED. USE RCOND IN CGECO FOR A RELIABLE
  39. C INDICATION OF SINGULARITY.
  40. C
  41. C LINPACK. THIS VERSION DATED 08/14/78 .
  42. C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  43. C
  44. C SUBROUTINES AND FUNCTIONS
  45. C
  46. C BLAS CAXPY,CSCAL,ICAMAX
  47. C FORTRAN ABS,AIMAG,REAL
  48. C
  49. C INTERNAL VARIABLES
  50. C
  51. COMPLEX T
  52. INTEGER ICAMAX,J,K,KP1,L,NM1
  53. C
  54. COMPLEX ZDUM
  55. REAL CABS1
  56. CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
  57. C
  58. C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
  59. C
  60. INFO = 0
  61. NM1 = N - 1
  62. IF (NM1 .LT. 1) GO TO 70
  63. DO 60 K = 1, NM1
  64. KP1 = K + 1
  65. C
  66. C FIND L = PIVOT INDEX
  67. C
  68. L = ICAMAX(N-K+1,A(K,K),1) + K - 1
  69. IPVT(K) = L
  70. C
  71. C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
  72. C
  73. IF (CABS1(A(L,K)) .EQ. 0.0E0) GO TO 40
  74. C
  75. C INTERCHANGE IF NECESSARY
  76. C
  77. IF (L .EQ. K) GO TO 10
  78. T = A(L,K)
  79. A(L,K) = A(K,K)
  80. A(K,K) = T
  81. 10 CONTINUE
  82. C
  83. C COMPUTE MULTIPLIERS
  84. C
  85. T = -(1.0E0,0.0E0)/A(K,K)
  86. CALL CSCAL(N-K,T,A(K+1,K),1)
  87. C
  88. C ROW ELIMINATION WITH COLUMN INDEXING
  89. C
  90. DO 30 J = KP1, N
  91. T = A(L,J)
  92. IF (L .EQ. K) GO TO 20
  93. A(L,J) = A(K,J)
  94. A(K,J) = T
  95. 20 CONTINUE
  96. CALL CAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
  97. 30 CONTINUE
  98. GO TO 50
  99. 40 CONTINUE
  100. INFO = K
  101. 50 CONTINUE
  102. 60 CONTINUE
  103. 70 CONTINUE
  104. IPVT(N) = N
  105. IF (CABS1(A(N,N)) .EQ. 0.0E0) INFO = N
  106. RETURN
  107. END