Hi all,
I wonder how can I form QR matrix after using qrdcmp() function.

The books says " In most applications we don't need to form Q explicitly"

Anybody knows how to do it by using NR? I am very frustrated, 'cause most of other packages using LAPACK form which is different from NR.

Many thanks.

Fang

If by that you mean you want to form the orthogonal matrix Q (R, the upper-triangular matrix, overwrites the upper part of the original matrix input in the qrdcmp() routine), you can use the other stuff output by the routine to form matrix Q, if you want it.

You must first realize that Q is not stored explicitly, but rather as a set of Householder vectors in the lower portion of the original matrix. Associated with that is the numbers in the vector c also output by the routine. (Read the book carefully for more details, especially chapter 11.)

It turns out that that accumulating these Householder transformations is quite easy! The answer to your question is implicit in the accompanying routine qrsolv()! Specifically, the first loop in the routine, applied to the columns of the identity matrix of corresponding dimension, will give the columns of matrix Q.

Hope this helps...

Jan M. (^_^)

P.S. The comment in the book was correct. Since storing Q as a sequence of Householder vectors is more compact, and that there is an easy way of applying the Householder transformations, you really would not need to form Q most of the time...

Thanks.

I understand the book and know want to save storage and computation, but just wonder if anyone writes one.

Anyway, I found a qr decomposiotion in GSL - The GNU Scientific Library. I copied their routine, problem solved.

I'd like to request this function in the next edition.

Fang
Best regards

ps. I also used the routine qr in matlab and I like the output implementation.

Hi

I have a problem recovering Q matrix explicitly from qrdcmp.

I have read carefully the explanations but something goes wrong
for instance with input matrix
ain = -0.9601171017 0.2795982659 0.2795982659 0.9601171017
which is already orthogonal,
I get for Q

-0.9601172209 1.279598236 1.279598236 0.9601171017

that is to say 1.27... instead of .27...

the matrix R seems correct
1.000000000 0.0000000000E+00 -0.2980232239E-07 1.000000000


here is what I do :
trying to compute the product Qdim-1....Q1


dim = size(a)
call qrdmcp(a,c,d,sing)
QR = Identity(dim) ! initialize for subsequent loop

do j = 1, dim-1
uj = 0.
uj(j: ) = a(j:,j)
QR = matmul(1 - &
(spread(uj, dim=2, ncopies = dim) * spread(uj, dim=1, ncopies = dim))/c(j),QR) ! outerproduct
end do


any help greatly appreciated
thanks

PS

the loop should be reversed to j=dim-1, 1,-1 for the correct Q1..Qn-1
matrix multiplication sequence, but anyway for dim=2
the output example is correct

finally it works the problem was I thought

that p1039, in
1 -u_j u_j' /c_j

I thought 1 was to be interpreted as a matrix with only ones, but it is the identity matrix