1 subroutine zchdc(a,lda,p,work,jpvt,job,info) 2 integer lda,p,jpvt(1),job,info 3 complex*16 a(lda,1),work(1) 4c 5c zchdc computes the cholesky decomposition of a positive definite 6c matrix. a pivoting option allows the user to estimate the 7c condition of a positive definite matrix or determine the rank 8c of a positive semidefinite matrix. 9c 10c on entry 11c 12c a complex*16(lda,p). 13c a contains the matrix whose decomposition is to 14c be computed. onlt the upper half of a need be stored. 15c the lower part of the array a is not referenced. 16c 17c lda integer. 18c lda is the leading dimension of the array a. 19c 20c p integer. 21c p is the order of the matrix. 22c 23c work complex*16. 24c work is a work array. 25c 26c jpvt integer(p). 27c jpvt contains integers that control the selection 28c of the pivot elements, if pivoting has been requested. 29c each diagonal element a(k,k) 30c is placed in one of three classes according to the 31c value of jpvt(k). 32c 33c if jpvt(k) .gt. 0, then x(k) is an initial 34c element. 35c 36c if jpvt(k) .eq. 0, then x(k) is a free element. 37c 38c if jpvt(k) .lt. 0, then x(k) is a final element. 39c 40c before the decomposition is computed, initial elements 41c are moved by symmetric row and column interchanges to 42c the beginning of the array a and final 43c elements to the end. both initial and final elements 44c are frozen in place during the computation and only 45c free elements are moved. at the k-th stage of the 46c reduction, if a(k,k) is occupied by a free element 47c it is interchanged with the largest free element 48c a(l,l) with l .ge. k. jpvt is not referenced if 49c job .eq. 0. 50c 51c job integer. 52c job is an integer that initiates column pivoting. 53c if job .eq. 0, no pivoting is done. 54c if job .ne. 0, pivoting is done. 55c 56c on return 57c 58c a a contains in its upper half the cholesky factor 59c of the matrix a as it has been permuted by pivoting. 60c 61c jpvt jpvt(j) contains the index of the diagonal element 62c of a that was moved into the j-th position, 63c provided pivoting was requested. 64c 65c info contains the index of the last positive diagonal 66c element of the cholesky factor. 67c 68c for positive definite matrices info = p is the normal return. 69c for pivoting with positive semidefinite matrices info will 70c in general be less than p. however, info may be greater than 71c the rank of a, since rounding error can cause an otherwise zero 72c element to be positive. indefinite systems will always cause 73c info to be less than p. 74c 75c linpack. this version dated 03/19/79 . 76c j.j. dongarra and g.w. stewart, argonne national laboratory and 77c university of maryland. 78c 79c 80c blas zaxpy,zswap 81c fortran dsqrt,dconjg 82c 83c internal variables 84c 85 integer pu,pl,plp1,i,j,jp,jt,k,kb,km1,kp1,l,maxl 86 complex*16 temp 87 double precision maxdia 88 logical swapk,negk 89 double precision dreal 90 complex*16 zdumr 91 dreal(zdumr) = zdumr 92c 93 pl = 1 94 pu = 0 95 info = p 96 if (job .eq. 0) go to 160 97c 98c pivoting has been requested. rearrange the 99c the elements according to jpvt. 100c 101 do 70 k = 1, p 102 swapk = jpvt(k) .gt. 0 103 negk = jpvt(k) .lt. 0 104 jpvt(k) = k 105 if (negk) jpvt(k) = -jpvt(k) 106 if (.not.swapk) go to 60 107 if (k .eq. pl) go to 50 108 call zswap(pl-1,a(1,k),1,a(1,pl),1) 109 temp = a(k,k) 110 a(k,k) = a(pl,pl) 111 a(pl,pl) = temp 112 a(pl,k) = dconjg(a(pl,k)) 113 plp1 = pl + 1 114 if (p .lt. plp1) go to 40 115 do 30 j = plp1, p 116 if (j .ge. k) go to 10 117 temp = dconjg(a(pl,j)) 118 a(pl,j) = dconjg(a(j,k)) 119 a(j,k) = temp 120 go to 20 121 10 continue 122 if (j .eq. k) go to 20 123 temp = a(k,j) 124 a(k,j) = a(pl,j) 125 a(pl,j) = temp 126 20 continue 127 30 continue 128 40 continue 129 jpvt(k) = jpvt(pl) 130 jpvt(pl) = k 131 50 continue 132 pl = pl + 1 133 60 continue 134 70 continue 135 pu = p 136 if (p .lt. pl) go to 150 137 do 140 kb = pl, p 138 k = p - kb + pl 139 if (jpvt(k) .ge. 0) go to 130 140 jpvt(k) = -jpvt(k) 141 if (pu .eq. k) go to 120 142 call zswap(k-1,a(1,k),1,a(1,pu),1) 143 temp = a(k,k) 144 a(k,k) = a(pu,pu) 145 a(pu,pu) = temp 146 a(k,pu) = dconjg(a(k,pu)) 147 kp1 = k + 1 148 if (p .lt. kp1) go to 110 149 do 100 j = kp1, p 150 if (j .ge. pu) go to 80 151 temp = dconjg(a(k,j)) 152 a(k,j) = dconjg(a(j,pu)) 153 a(j,pu) = temp 154 go to 90 155 80 continue 156 if (j .eq. pu) go to 90 157 temp = a(k,j) 158 a(k,j) = a(pu,j) 159 a(pu,j) = temp 160 90 continue 161 100 continue 162 110 continue 163 jt = jpvt(k) 164 jpvt(k) = jpvt(pu) 165 jpvt(pu) = jt 166 120 continue 167 pu = pu - 1 168 130 continue 169 140 continue 170 150 continue 171 160 continue 172 do 270 k = 1, p 173c 174c reduction loop. 175c 176 maxdia = dreal(a(k,k)) 177 kp1 = k + 1 178 maxl = k 179c 180c determine the pivot element. 181c 182 if (k .lt. pl .or. k .ge. pu) go to 190 183 do 180 l = kp1, pu 184 if (dreal(a(l,l)) .le. maxdia) go to 170 185 maxdia = dreal(a(l,l)) 186 maxl = l 187 170 continue 188 180 continue 189 190 continue 190c 191c quit if the pivot element is not positive. 192c 193 if (maxdia .gt. 0.0d0) go to 200 194 info = k - 1 195c ......exit 196 go to 280 197 200 continue 198 if (k .eq. maxl) go to 210 199c 200c start the pivoting and update jpvt. 201c 202 km1 = k - 1 203 call zswap(km1,a(1,k),1,a(1,maxl),1) 204 a(maxl,maxl) = a(k,k) 205 a(k,k) = dcmplx(maxdia,0.0d0) 206 jp = jpvt(maxl) 207 jpvt(maxl) = jpvt(k) 208 jpvt(k) = jp 209 a(k,maxl) = dconjg(a(k,maxl)) 210 210 continue 211c 212c reduction step. pivoting is contained across the rows. 213c 214 work(k) = dcmplx(dsqrt(dreal(a(k,k))),0.0d0) 215 a(k,k) = work(k) 216 if (p .lt. kp1) go to 260 217 do 250 j = kp1, p 218 if (k .eq. maxl) go to 240 219 if (j .ge. maxl) go to 220 220 temp = dconjg(a(k,j)) 221 a(k,j) = dconjg(a(j,maxl)) 222 a(j,maxl) = temp 223 go to 230 224 220 continue 225 if (j .eq. maxl) go to 230 226 temp = a(k,j) 227 a(k,j) = a(maxl,j) 228 a(maxl,j) = temp 229 230 continue 230 240 continue 231 a(k,j) = a(k,j)/work(k) 232 work(j) = dconjg(a(k,j)) 233 temp = -a(k,j) 234 call zaxpy(j-k,temp,work(kp1),1,a(kp1,j),1) 235 250 continue 236 260 continue 237 270 continue 238 280 continue 239 return 240 end 241