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