subroutine cqrdc(x,ldx,n,p,qraux,jpvt,work,job)
integer ldx,n,p,job
integer jpvt(1)
complex x(ldx,1),qraux(1),work(1)
c
c cqrdc uses householder transformations to compute the qr
c factorization of an n by p matrix x. column pivoting
c based on the 2-norms of the reduced columns may be
c performed at the users option.
c
c on entry
c
c x complex(ldx,p), where ldx .ge. n.
c x contains the matrix whose decomposition is to be
c computed.
c
c ldx integer.
c ldx is the leading dimension of the array x.
c
c n integer.
c n is the number of rows of the matrix x.
c
c p integer.
c p is the number of columns of the matrix x.
c
c jpvt integer(p).
c jpvt contains integers that control the selection
c of the pivot columns. the k-th column x(k) of x
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 column.
c
c if jpvt(k) .eq. 0, then x(k) is a free column.
c
c if jpvt(k) .lt. 0, then x(k) is a final column.
c
c before the decomposition is computed, initial columns
c are moved to the beginning of the array x and final
c columns to the end. both initial and final columns
c are frozen in place during the computation and only
c free columns are moved. at the k-th stage of the
c reduction, if x(k) is occupied by a free column
c it is interchanged with the free column of largest
c reduced norm. jpvt is not referenced if
c job .eq. 0.
c
c work complex(p).
c work is a work array. work 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 x x contains in its upper triangle the upper
c triangular matrix r of the qr factorization.
c below its diagonal x contains information from
c which the unitary part of the decomposition
c can be recovered. note that if pivoting has
c been requested, the decomposition is not that
c of the original matrix x but that of x
c with its columns permuted as described by jpvt.
c
c qraux complex(p).
c qraux contains further information required to recover
c the unitary part of the decomposition.
c
c jpvt jpvt(k) contains the index of the column of the
c original matrix that has been interchanged into
c the k-th column, if pivoting was requested.
c
c linpack. this version dated 08/14/78 .
c g.w. stewart, university of maryland, argonne national lab.
c
c cqrdc uses the following functions and subprograms.
c
c blas caxpy,cdotc,cscal,cswap,scnrm2
c fortran abs,aimag,amax1,cabs,cmplx,csqrt,min0,real
c
c internal variables
c
integer j,jp,l,lp1,lup,maxj,pl,pu
real maxnrm,scnrm2,tt
complex cdotc,nrmxl,t
logical negj,swapj
c
complex csign,zdum,zdum1,zdum2
real cabs1
csign(zdum1,zdum2) = cabs(zdum1)*(zdum2/cabs(zdum2))
cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum))
c
pl = 1
pu = 0
if (job .eq. 0) go to 60
c
c pivoting has been requested. rearrange the columns
c according to jpvt.
c
do 20 j = 1, p
swapj = jpvt(j) .gt. 0
negj = jpvt(j) .lt. 0
jpvt(j) = j
if (negj) jpvt(j) = -j
if (.not.swapj) go to 10
if (j .ne. pl) call cswap(n,x(1,pl),1,x(1,j),1)
jpvt(j) = jpvt(pl)
jpvt(pl) = j
pl = pl + 1
10 continue
20 continue
pu = p
do 50 jj = 1, p
j = p - jj + 1
if (jpvt(j) .ge. 0) go to 40
jpvt(j) = -jpvt(j)
if (j .eq. pu) go to 30
call cswap(n,x(1,pu),1,x(1,j),1)
jp = jpvt(pu)
jpvt(pu) = jpvt(j)
jpvt(j) = jp
30 continue
pu = pu - 1
40 continue
50 continue
60 continue
c
c compute the norms of the free columns.
c
if (pu .lt. pl) go to 80
do 70 j = pl, pu
qraux(j) = cmplx(scnrm2(n,x(1,j),1),0.0e0)
work(j) = qraux(j)
70 continue
80 continue
c
c perform the householder reduction of x.
c
lup = min0(n,p)
do 200 l = 1, lup
if (l .lt. pl .or. l .ge. pu) go to 120
c
c locate the column of largest norm and bring it
c into the pivot position.
c
maxnrm = 0.0e0
maxj = l
do 100 j = l, pu
if (real(qraux(j)) .le. maxnrm) go to 90
maxnrm = real(qraux(j))
maxj = j
90 continue
100 continue
if (maxj .eq. l) go to 110
call cswap(n,x(1,l),1,x(1,maxj),1)
qraux(maxj) = qraux(l)
work(maxj) = work(l)
jp = jpvt(maxj)
jpvt(maxj) = jpvt(l)
jpvt(l) = jp
110 continue
120 continue
qraux(l) = (0.0e0,0.0e0)
if (l .eq. n) go to 190
c
c compute the householder transformation for column l.
c
nrmxl = cmplx(scnrm2(n-l+1,x(l,l),1),0.0e0)
if (cabs1(nrmxl) .eq. 0.0e0) go to 180
if (cabs1(x(l,l)) .ne. 0.0e0)
* nrmxl = csign(nrmxl,x(l,l))
call cscal(n-l+1,(1.0e0,0.0e0)/nrmxl,x(l,l),1)
x(l,l) = (1.0e0,0.0e0) + x(l,l)
c
c apply the transformation to the remaining columns,
c updating the norms.
c
lp1 = l + 1
if (p .lt. lp1) go to 170
do 160 j = lp1, p
t = -cdotc(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
call caxpy(n-l+1,t,x(l,l),1,x(l,j),1)
if (j .lt. pl .or. j .gt. pu) go to 150
if (cabs1(qraux(j)) .eq. 0.0e0) go to 150
tt = 1.0e0 - (cabs(x(l,j))/real(qraux(j)))**2
tt = amax1(tt,0.0e0)
t = cmplx(tt,0.0e0)
tt = 1.0e0
* + 0.05e0*tt*(real(qraux(j))/real(work(j)))**2
if (tt .eq. 1.0e0) go to 130
qraux(j) = qraux(j)*csqrt(t)
go to 140
130 continue
qraux(j) = cmplx(scnrm2(n-l,x(l+1,j),1),0.0e0)
work(j) = qraux(j)
140 continue
150 continue
160 continue
170 continue
c
c save the transformation.
c
qraux(l) = x(l,l)
x(l,l) = -nrmxl
180 continue
190 continue
200 continue
return
end