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
.