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