subroutine cgefa(a,lda,n,ipvt,info) integer lda,n,ipvt(1),info complex a(lda,1) c c cgefa factors a complex matrix by gaussian elimination. c c cgefa is usually called by cgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for cgeco) = (1 + 9/n)*(time for cgefa) . c c on entry c c a complex(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that cgesl or cgedi will divide by zero c if called. use rcond in cgeco for a reliable c indication of singularity. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas caxpy,cscal,icamax c fortran abs,aimag,real c c internal variables c complex t integer icamax,j,k,kp1,l,nm1 c complex zdum real cabs1 cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum)) c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = icamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (cabs1(a(l,k)) .eq. 0.0e0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -(1.0e0,0.0e0)/a(k,k) call cscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call caxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (cabs1(a(n,n)) .eq. 0.0e0) info = n return end