[SciPy-dev] RQ Decomposition
Harry Kalogirou
harkal at sylphis3d.com
Fri Jul 22 04:02:28 CDT 2005
Hello!
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 :
<----code----code---->
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)
<----code----code---->
It does not work. As I don't understand fortran indexing scheme
maybe i'm doing something wrong...
can anyone help??
Best regards,
Harry
More information about the Scipy-dev
mailing list