subroutine cqrz(a,q,s,m,l,ldq,lds)
      integer m,l,ldq,lds
      double complex a(m),q(ldq,l),s(lds,l)
c
c     cqrz computes the qr factorization in the form
c     a * r(inverse) = q
c     of the double complex column-circulant matrix a .
c
c     on entry
c
c        a       double complex(m)
c                the first column of the column-circulant matrix .
c                on return a is unaltered .
c
c        m       integer
c                the number of rows of the matrices a and q .
c                m must be at least as large as l .
c
c        l       integer
c                the number of columns of the matrices a and q
c                and the order of the upper triangular matrix s .
c
c        ldq     integer
c                the leading dimension of the array q .
c
c        lds     integer
c                the leading dimension of the array s .
c
c     on return
c
c        q       double complex(m,l)
c                the q matrix of the factorization .
c                the columns of q are orthonormal .
c
c        s       double complex(l,l)
c                the inverse of the r matrix of the factorization .
c                elements below the main diagonal are not accessed .
c
c     toeplitz package. this version dated 07/23/82 .
c
c     subroutines and functions
c
c        linpack ... zaxpy,zdotc,zdscal,dznrm2
c        fortran ... dconjg
c
c     internal variables
c
      integer i,j,j1,ji
      double precision scale,dznrm2
      double complex c,zdotc
c
c     initialization (last column of q used as work vector) .
c
      do 10 i = 1, m
         q(i,1) = a(i)
         q(i,l) = a(i)
   10 continue
c
c     recurrent process for the lattice algorithm with normalization .
c
      do 70 j1 = 1, l
         j = j1 + 1
         scale = 1.0d0/dznrm2(m,q(1,j1),1)
         if (j1 .eq. l) go to 60
            c = -scale*(dconjg(q(m,j1))*q(1,l) +
     *           zdotc(m-1,q(1,j1),1,q(2,l),1))/dznrm2(m,q(1,l),1)
            q(1,j) = q(m,j1) + c*q(1,l)
            do 20 i = 2, m
               q(i,j) = q(i-1,j1) + c*q(i,l)
   20       continue
            if (j .eq. l) go to 30
               q(1,l) = q(1,l) + c*q(m,j1)
               call zaxpy(m-1,c,q(1,j1),1,q(2,l),1)
   30       continue
            s(1,j) = c
            if (j .eq. 2) go to 50
               do 40 i = 2, j1
                  ji = j - i
                  s(i,j) = s(i-1,j1) + c*s(ji,j1)
   40          continue
   50       continue
   60    continue
         call zdscal(m,scale,q(1,j1),1)
         s(j1,j1) = (1.0d0,0.0d0)
         call zdscal(j1,scale,s(1,j1),1)
   70 continue
      return
      end