[SciPy-dev] RQ Decomposition

Harry Kalogirou harkal at sylphis3d.com
Fri Jul 22 04:02:28 CDT 2005


I'm trying to implement th RQ decomposition in SciPy. Unfortunatly 
I'm not good at fortran and I have some problems.

So far I created the wrapper function for "dgerqf" based on the
already provided wrapper for the "dgeqrf" function.

Now I try to implement the "linalg.rq" function based on the
"linalg.qr". The problem is that I can't compute Q as said in the 
comments of the "dgerqf" function. 

>From the comments of the "dgerqf" function :

>  where tau is a real scalar, and v is a real vector with
>  v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in
>  A(m-k+i,1:n-k+i-1), and tau in TAU(i).

The code I write is :

    k = min(M,N)
    for i in range(k):
        v = zeros((M,),t)
        v[N-k+i] = 1
        v[0:N-k+i-1] = rq[M-k+i,0:N-k+i-1]
        H = gemm(-tau[i],v,v,1+0j,ident,trans_b=2)
        Q = gemm(1,Q,H)

It does not work. As I don't understand fortran indexing scheme
maybe i'm doing something wrong...

can anyone help??

Best regards,

More information about the Scipy-dev mailing list