1 subroutine chifa(a,lda,n,kpvt,info) 2 integer lda,n,kpvt(1),info 3 complex a(lda,1) 4c 5c chifa factors a complex hermitian matrix by elimination 6c with symmetric pivoting. 7c 8c to solve a*x = b , follow chifa by chisl. 9c to compute inverse(a)*c , follow chifa by chisl. 10c to compute determinant(a) , follow chifa by chidi. 11c to compute inertia(a) , follow chifa by chidi. 12c to compute inverse(a) , follow chifa by chidi. 13c 14c on entry 15c 16c a complex(lda,n) 17c the hermitian matrix to be factored. 18c only the diagonal and upper triangle are used. 19c 20c lda integer 21c the leading dimension of the array a . 22c 23c n integer 24c the order of the matrix a . 25c 26c on return 27c 28c a a block diagonal matrix and the multipliers which 29c were used to obtain it. 30c the factorization can be written a = u*d*ctrans(u) 31c where u is a product of permutation and unit 32c upper triangular matrices , ctrans(u) is the 33c conjugate transpose of u , and d is block diagonal 34c with 1 by 1 and 2 by 2 blocks. 35c 36c kpvt integer(n) 37c an integer vector of pivot indices. 38c 39c info integer 40c = 0 normal value. 41c = k if the k-th pivot block is singular. this is 42c not an error condition for this subroutine, 43c but it does indicate that chisl or chidi may 44c divide by zero if called. 45c 46c linpack. this version dated 08/14/78 . 47c james bunch, univ. calif. san diego, argonne nat. lab. 48c 49c subroutines and functions 50c 51c blas caxpy,cswap,icamax 52c fortran abs,aimag,amax1,cmplx,conjg,real,sqrt 53c 54c internal variables 55c 56 complex ak,akm1,bk,bkm1,denom,mulk,mulkm1,t 57 real absakk,alpha,colmax,rowmax 58 integer imax,imaxp1,j,jj,jmax,k,km1,km2,kstep,icamax 59 logical swap 60c 61 complex zdum 62 real cabs1 63 cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum)) 64c 65c initialize 66c 67c alpha is used in choosing pivot block size. 68 alpha = (1.0e0 + sqrt(17.0e0))/8.0e0 69c 70 info = 0 71c 72c main loop on k, which goes from n to 1. 73c 74 k = n 75 10 continue 76c 77c leave the loop if k=0 or k=1. 78c 79c ...exit 80 if (k .eq. 0) go to 200 81 if (k .gt. 1) go to 20 82 kpvt(1) = 1 83 if (cabs1(a(1,1)) .eq. 0.0e0) info = 1 84c ......exit 85 go to 200 86 20 continue 87c 88c this section of code determines the kind of 89c elimination to be performed. when it is completed, 90c kstep will be set to the size of the pivot block, and 91c swap will be set to .true. if an interchange is 92c required. 93c 94 km1 = k - 1 95 absakk = cabs1(a(k,k)) 96c 97c determine the largest off-diagonal element in 98c column k. 99c 100 imax = icamax(k-1,a(1,k),1) 101 colmax = cabs1(a(imax,k)) 102 if (absakk .lt. alpha*colmax) go to 30 103 kstep = 1 104 swap = .false. 105 go to 90 106 30 continue 107c 108c determine the largest off-diagonal element in 109c row imax. 110c 111 rowmax = 0.0e0 112 imaxp1 = imax + 1 113 do 40 j = imaxp1, k 114 rowmax = amax1(rowmax,cabs1(a(imax,j))) 115 40 continue 116 if (imax .eq. 1) go to 50 117 jmax = icamax(imax-1,a(1,imax),1) 118 rowmax = amax1(rowmax,cabs1(a(jmax,imax))) 119 50 continue 120 if (cabs1(a(imax,imax)) .lt. alpha*rowmax) go to 60 121 kstep = 1 122 swap = .true. 123 go to 80 124 60 continue 125 if (absakk .lt. alpha*colmax*(colmax/rowmax)) go to 70 126 kstep = 1 127 swap = .false. 128 go to 80 129 70 continue 130 kstep = 2 131 swap = imax .ne. km1 132 80 continue 133 90 continue 134 if (amax1(absakk,colmax) .ne. 0.0e0) go to 100 135c 136c column k is zero. set info and iterate the loop. 137c 138 kpvt(k) = k 139 info = k 140 go to 190 141 100 continue 142 if (kstep .eq. 2) go to 140 143c 144c 1 x 1 pivot block. 145c 146 if (.not.swap) go to 120 147c 148c perform an interchange. 149c 150 call cswap(imax,a(1,imax),1,a(1,k),1) 151 do 110 jj = imax, k 152 j = k + imax - jj 153 t = conjg(a(j,k)) 154 a(j,k) = conjg(a(imax,j)) 155 a(imax,j) = t 156 110 continue 157 120 continue 158c 159c perform the elimination. 160c 161 do 130 jj = 1, km1 162 j = k - jj 163 mulk = -a(j,k)/a(k,k) 164 t = conjg(mulk) 165 call caxpy(j,t,a(1,k),1,a(1,j),1) 166 a(j,j) = cmplx(real(a(j,j)),0.0e0) 167 a(j,k) = mulk 168 130 continue 169c 170c set the pivot array. 171c 172 kpvt(k) = k 173 if (swap) kpvt(k) = imax 174 go to 190 175 140 continue 176c 177c 2 x 2 pivot block. 178c 179 if (.not.swap) go to 160 180c 181c perform an interchange. 182c 183 call cswap(imax,a(1,imax),1,a(1,k-1),1) 184 do 150 jj = imax, km1 185 j = km1 + imax - jj 186 t = conjg(a(j,k-1)) 187 a(j,k-1) = conjg(a(imax,j)) 188 a(imax,j) = t 189 150 continue 190 t = a(k-1,k) 191 a(k-1,k) = a(imax,k) 192 a(imax,k) = t 193 160 continue 194c 195c perform the elimination. 196c 197 km2 = k - 2 198 if (km2 .eq. 0) go to 180 199 ak = a(k,k)/a(k-1,k) 200 akm1 = a(k-1,k-1)/conjg(a(k-1,k)) 201 denom = 1.0e0 - ak*akm1 202 do 170 jj = 1, km2 203 j = km1 - jj 204 bk = a(j,k)/a(k-1,k) 205 bkm1 = a(j,k-1)/conjg(a(k-1,k)) 206 mulk = (akm1*bk - bkm1)/denom 207 mulkm1 = (ak*bkm1 - bk)/denom 208 t = conjg(mulk) 209 call caxpy(j,t,a(1,k),1,a(1,j),1) 210 t = conjg(mulkm1) 211 call caxpy(j,t,a(1,k-1),1,a(1,j),1) 212 a(j,k) = mulk 213 a(j,k-1) = mulkm1 214 a(j,j) = cmplx(real(a(j,j)),0.0e0) 215 170 continue 216 180 continue 217c 218c set the pivot array. 219c 220 kpvt(k) = 1 - k 221 if (swap) kpvt(k) = -imax 222 kpvt(k-1) = kpvt(k) 223 190 continue 224 k = k - kstep 225 go to 10 226 200 continue 227 return 228 end 229