subroutine cqrs(a,q,s,m,l,ldq,lds) integer m,l,ldq,lds real a(m),q(ldq,l),s(lds,l) c c cqrs computes the qr factorization in the form c a * r(inverse) = q c of the real column-circulant matrix a . c c on entry c c a real(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 real(m,l) c the q matrix of the factorization . c the columns of q are orthonormal . c c s real(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 ... saxpy,sdot,sscal,snrm2 c c internal variables c integer i,j,j1,ji real scale,snrm2 real c,sdot 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.0e0/snrm2(m,q(1,j1),1) if (j1 .eq. l) go to 60 c = -scale*(q(m,j1)*q(1,l) + * sdot(m-1,q(1,j1),1,q(2,l),1))/snrm2(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 saxpy(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 sscal(m,scale,q(1,j1),1) s(j1,j1) = 1.0e0 call sscal(j1,scale,s(1,j1),1) 70 continue return end