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