1      subroutine cqrdc(x,ldx,n,p,qraux,jpvt,work,job)
2      integer ldx,n,p,job
3      integer jpvt(1)
4      complex x(ldx,1),qraux(1),work(1)
5c
6c     cqrdc uses householder transformations to compute the qr
7c     factorization of an n by p matrix x.  column pivoting
8c     based on the 2-norms of the reduced columns may be
9c     performed at the users option.
10c
11c     on entry
12c
13c        x       complex(ldx,p), where ldx .ge. n.
14c                x contains the matrix whose decomposition is to be
15c                computed.
16c
17c        ldx     integer.
18c                ldx is the leading dimension of the array x.
19c
20c        n       integer.
21c                n is the number of rows of the matrix x.
22c
23c        p       integer.
24c                p is the number of columns of the matrix x.
25c
26c        jpvt    integer(p).
27c                jpvt contains integers that control the selection
28c                of the pivot columns.  the k-th column x(k) of x
29c                is placed in one of three classes according to the
30c                value of jpvt(k).
31c
32c                   if jpvt(k) .gt. 0, then x(k) is an initial
33c                                      column.
34c
35c                   if jpvt(k) .eq. 0, then x(k) is a free column.
36c
37c                   if jpvt(k) .lt. 0, then x(k) is a final column.
38c
39c                before the decomposition is computed, initial columns
40c                are moved to the beginning of the array x and final
41c                columns to the end.  both initial and final columns
42c                are frozen in place during the computation and only
43c                free columns are moved.  at the k-th stage of the
44c                reduction, if x(k) is occupied by a free column
45c                it is interchanged with the free column of largest
46c                reduced norm.  jpvt is not referenced if
47c                job .eq. 0.
48c
49c        work    complex(p).
50c                work is a work array.  work is not referenced if
51c                job .eq. 0.
52c
53c        job     integer.
54c                job is an integer that initiates column pivoting.
55c                if job .eq. 0, no pivoting is done.
56c                if job .ne. 0, pivoting is done.
57c
58c     on return
59c
60c        x       x contains in its upper triangle the upper
61c                triangular matrix r of the qr factorization.
62c                below its diagonal x contains information from
63c                which the unitary part of the decomposition
64c                can be recovered.  note that if pivoting has
65c                been requested, the decomposition is not that
66c                of the original matrix x but that of x
67c                with its columns permuted as described by jpvt.
68c
69c        qraux   complex(p).
70c                qraux contains further information required to recover
71c                the unitary part of the decomposition.
72c
73c        jpvt    jpvt(k) contains the index of the column of the
74c                original matrix that has been interchanged into
75c                the k-th column, if pivoting was requested.
76c
77c     linpack. this version dated 08/14/78 .
78c     g.w. stewart, university of maryland, argonne national lab.
79c
80c     cqrdc uses the following functions and subprograms.
81c
82c     blas caxpy,cdotc,cscal,cswap,scnrm2
83c     fortran abs,aimag,amax1,cabs,cmplx,csqrt,min0,real
84c
85c     internal variables
86c
87      integer j,jp,l,lp1,lup,maxj,pl,pu
88      real maxnrm,scnrm2,tt
89      complex cdotc,nrmxl,t
90      logical negj,swapj
91c
92      complex csign,zdum,zdum1,zdum2
93      real cabs1
94      csign(zdum1,zdum2) = cabs(zdum1)*(zdum2/cabs(zdum2))
95      cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum))
96c
97      pl = 1
98      pu = 0
99      if (job .eq. 0) go to 60
100c
101c        pivoting has been requested.  rearrange the columns
102c        according to jpvt.
103c
104         do 20 j = 1, p
105            swapj = jpvt(j) .gt. 0
106            negj = jpvt(j) .lt. 0
107            jpvt(j) = j
108            if (negj) jpvt(j) = -j
109            if (.not.swapj) go to 10
110               if (j .ne. pl) call cswap(n,x(1,pl),1,x(1,j),1)
111               jpvt(j) = jpvt(pl)
112               jpvt(pl) = j
113               pl = pl + 1
114   10       continue
115   20    continue
116         pu = p
117         do 50 jj = 1, p
118            j = p - jj + 1
119            if (jpvt(j) .ge. 0) go to 40
120               jpvt(j) = -jpvt(j)
121               if (j .eq. pu) go to 30
122                  call cswap(n,x(1,pu),1,x(1,j),1)
123                  jp = jpvt(pu)
124                  jpvt(pu) = jpvt(j)
125                  jpvt(j) = jp
126   30          continue
127               pu = pu - 1
128   40       continue
129   50    continue
130   60 continue
131c
132c     compute the norms of the free columns.
133c
134      if (pu .lt. pl) go to 80
135      do 70 j = pl, pu
136         qraux(j) = cmplx(scnrm2(n,x(1,j),1),0.0e0)
137         work(j) = qraux(j)
138   70 continue
139   80 continue
140c
141c     perform the householder reduction of x.
142c
143      lup = min0(n,p)
144      do 200 l = 1, lup
145         if (l .lt. pl .or. l .ge. pu) go to 120
146c
147c           locate the column of largest norm and bring it
148c           into the pivot position.
149c
150            maxnrm = 0.0e0
151            maxj = l
152            do 100 j = l, pu
153               if (real(qraux(j)) .le. maxnrm) go to 90
154                  maxnrm = real(qraux(j))
155                  maxj = j
156   90          continue
157  100       continue
158            if (maxj .eq. l) go to 110
159               call cswap(n,x(1,l),1,x(1,maxj),1)
160               qraux(maxj) = qraux(l)
161               work(maxj) = work(l)
162               jp = jpvt(maxj)
163               jpvt(maxj) = jpvt(l)
164               jpvt(l) = jp
165  110       continue
166  120    continue
167         qraux(l) = (0.0e0,0.0e0)
168         if (l .eq. n) go to 190
169c
170c           compute the householder transformation for column l.
171c
172            nrmxl = cmplx(scnrm2(n-l+1,x(l,l),1),0.0e0)
173            if (cabs1(nrmxl) .eq. 0.0e0) go to 180
174               if (cabs1(x(l,l)) .ne. 0.0e0)
175     *            nrmxl = csign(nrmxl,x(l,l))
176               call cscal(n-l+1,(1.0e0,0.0e0)/nrmxl,x(l,l),1)
177               x(l,l) = (1.0e0,0.0e0) + x(l,l)
178c
179c              apply the transformation to the remaining columns,
180c              updating the norms.
181c
182               lp1 = l + 1
183               if (p .lt. lp1) go to 170
184               do 160 j = lp1, p
185                  t = -cdotc(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
186                  call caxpy(n-l+1,t,x(l,l),1,x(l,j),1)
187                  if (j .lt. pl .or. j .gt. pu) go to 150
188                  if (cabs1(qraux(j)) .eq. 0.0e0) go to 150
189                     tt = 1.0e0 - (cabs(x(l,j))/real(qraux(j)))**2
190                     tt = amax1(tt,0.0e0)
191                     t = cmplx(tt,0.0e0)
192                     tt = 1.0e0
193     *                    + 0.05e0*tt*(real(qraux(j))/real(work(j)))**2
194                     if (tt .eq. 1.0e0) go to 130
195                        qraux(j) = qraux(j)*csqrt(t)
196                     go to 140
197  130                continue
198                        qraux(j) = cmplx(scnrm2(n-l,x(l+1,j),1),0.0e0)
199                        work(j) = qraux(j)
200  140                continue
201  150             continue
202  160          continue
203  170          continue
204c
205c              save the transformation.
206c
207               qraux(l) = x(l,l)
208               x(l,l) = -nrmxl
209  180       continue
210  190    continue
211  200 continue
212      return
213      end
214