| subroutine dchdc(a,lda,p,work,jpvt,job,info) |
| integer lda,p,jpvt(p),job,info |
| double precision a(lda,p),work(p) |
| c |
| c dchdc computes the cholesky decomposition of a positive definite |
| c matrix. a pivoting option allows the user to estimate the |
| c condition of a positive definite matrix or determine the rank |
| c of a positive semidefinite matrix. |
| c |
| c on entry |
| c |
| c a double precision(lda,p). |
| c a contains the matrix whose decomposition is to |
| c be computed. onlt the upper half of a need be stored. |
| c the lower part of the array a is not referenced. |
| c |
| c lda integer. |
| c lda is the leading dimension of the array a. |
| c |
| c p integer. |
| c p is the order of the matrix. |
| c |
| c work double precision. |
| c work is a work array. |
| c |
| c jpvt integer(p). |
| c jpvt contains integers that control the selection |
| c of the pivot elements, if pivoting has been requested. |
| c each diagonal element a(k,k) |
| c is placed in one of three classes according to the |
| c value of jpvt(k). |
| c |
| c if jpvt(k) .gt. 0, then x(k) is an initial |
| c element. |
| c |
| c if jpvt(k) .eq. 0, then x(k) is a free element. |
| c |
| c if jpvt(k) .lt. 0, then x(k) is a final element. |
| c |
| c before the decomposition is computed, initial elements |
| c are moved by symmetric row and column interchanges to |
| c the beginning of the array a and final |
| c elements to the end. both initial and final elements |
| c are frozen in place during the computation and only |
| c free elements are moved. at the k-th stage of the |
| c reduction, if a(k,k) is occupied by a free element |
| c it is interchanged with the largest free element |
| c a(l,l) with l .ge. k. jpvt is not referenced if |
| c job .eq. 0. |
| c |
| c job integer. |
| c job is an integer that initiates column pivoting. |
| c if job .eq. 0, no pivoting is done. |
| c if job .ne. 0, pivoting is done. |
| c |
| c on return |
| c |
| c a a contains in its upper half the cholesky factor |
| c of the matrix a as it has been permuted by pivoting. |
| c |
| c jpvt jpvt(j) contains the index of the diagonal element |
| c of a that was moved into the j-th position, |
| c provided pivoting was requested. |
| c |
| c info contains the index of the last positive diagonal |
| c element of the cholesky factor. |
| c |
| c for positive definite matrices info = p is the normal return. |
| c for pivoting with positive semidefinite matrices info will |
| c in general be less than p. however, info may be greater than |
| c the rank of a, since rounding error can cause an otherwise zero |
| c element to be positive. indefinite systems will always cause |
| c info to be less than p. |
| c |
| c linpack. this version dated 08/14/78 . |
| c j.j. dongarra and g.w. stewart, argonne national laboratory and |
| c university of maryland. |
| c |
| c |
| c blas daxpy,dswap |
| c fortran dsqrt |
| c |
| c internal variables |
| c |
| integer pu,pl,plp1,j,jp,jt,k,kb,km1,kp1,l,maxl |
| double precision temp |
| double precision maxdia |
| logical swapk,negk |
| c |
| pl = 1 |
| pu = 0 |
| info = p |
| if (job .eq. 0) go to 160 |
| c |
| c pivoting has been requested. rearrange the |
| c the elements according to jpvt. |
| c |
| do 70 k = 1, p |
| swapk = jpvt(k) .gt. 0 |
| negk = jpvt(k) .lt. 0 |
| jpvt(k) = k |
| if (negk) jpvt(k) = -jpvt(k) |
| if (.not.swapk) go to 60 |
| if (k .eq. pl) go to 50 |
| call dswap(pl-1,a(1,k),1,a(1,pl),1) |
| temp = a(k,k) |
| a(k,k) = a(pl,pl) |
| a(pl,pl) = temp |
| plp1 = pl + 1 |
| if (p .lt. plp1) go to 40 |
| do 30 j = plp1, p |
| if (j .ge. k) go to 10 |
| temp = a(pl,j) |
| a(pl,j) = a(j,k) |
| a(j,k) = temp |
| go to 20 |
| 10 continue |
| if (j .eq. k) go to 20 |
| temp = a(k,j) |
| a(k,j) = a(pl,j) |
| a(pl,j) = temp |
| 20 continue |
| 30 continue |
| 40 continue |
| jpvt(k) = jpvt(pl) |
| jpvt(pl) = k |
| 50 continue |
| pl = pl + 1 |
| 60 continue |
| 70 continue |
| pu = p |
| if (p .lt. pl) go to 150 |
| do 140 kb = pl, p |
| k = p - kb + pl |
| if (jpvt(k) .ge. 0) go to 130 |
| jpvt(k) = -jpvt(k) |
| if (pu .eq. k) go to 120 |
| call dswap(k-1,a(1,k),1,a(1,pu),1) |
| temp = a(k,k) |
| a(k,k) = a(pu,pu) |
| a(pu,pu) = temp |
| kp1 = k + 1 |
| if (p .lt. kp1) go to 110 |
| do 100 j = kp1, p |
| if (j .ge. pu) go to 80 |
| temp = a(k,j) |
| a(k,j) = a(j,pu) |
| a(j,pu) = temp |
| go to 90 |
| 80 continue |
| if (j .eq. pu) go to 90 |
| temp = a(k,j) |
| a(k,j) = a(pu,j) |
| a(pu,j) = temp |
| 90 continue |
| 100 continue |
| 110 continue |
| jt = jpvt(k) |
| jpvt(k) = jpvt(pu) |
| jpvt(pu) = jt |
| 120 continue |
| pu = pu - 1 |
| 130 continue |
| 140 continue |
| 150 continue |
| 160 continue |
| do 270 k = 1, p |
| c |
| c reduction loop. |
| c |
| maxdia = a(k,k) |
| kp1 = k + 1 |
| maxl = k |
| c |
| c determine the pivot element. |
| c |
| if (k .lt. pl .or. k .ge. pu) go to 190 |
| do 180 l = kp1, pu |
| if (a(l,l) .le. maxdia) go to 170 |
| maxdia = a(l,l) |
| maxl = l |
| 170 continue |
| 180 continue |
| 190 continue |
| c |
| c quit if the pivot element is not positive. |
| c |
| if (maxdia .gt. 0.0d0) go to 200 |
| info = k - 1 |
| c ......exit |
| go to 280 |
| 200 continue |
| if (k .eq. maxl) go to 210 |
| c |
| c start the pivoting and update jpvt. |
| c |
| km1 = k - 1 |
| call dswap(km1,a(1,k),1,a(1,maxl),1) |
| a(maxl,maxl) = a(k,k) |
| a(k,k) = maxdia |
| jp = jpvt(maxl) |
| jpvt(maxl) = jpvt(k) |
| jpvt(k) = jp |
| 210 continue |
| c |
| c reduction step. pivoting is contained across the rows. |
| c |
| work(k) = dsqrt(a(k,k)) |
| a(k,k) = work(k) |
| if (p .lt. kp1) go to 260 |
| do 250 j = kp1, p |
| if (k .eq. maxl) go to 240 |
| if (j .ge. maxl) go to 220 |
| temp = a(k,j) |
| a(k,j) = a(j,maxl) |
| a(j,maxl) = temp |
| go to 230 |
| 220 continue |
| if (j .eq. maxl) go to 230 |
| temp = a(k,j) |
| a(k,j) = a(maxl,j) |
| a(maxl,j) = temp |
| 230 continue |
| 240 continue |
| a(k,j) = a(k,j)/work(k) |
| work(j) = a(k,j) |
| temp = -a(k,j) |
| call daxpy(j-k,temp,work(kp1),1,a(kp1,j),1) |
| 250 continue |
| 260 continue |
| 270 continue |
| 280 continue |
| return |
| end |
