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