0

I'm trying to get Q-Matrix from LAPACK zgeqrf and zungqr routines.

I have a Nw-by-Nw complex matrix with non-orthogonal vectors on its columns. This is my C++ code. (The matrix is named vr_tr)

//QR Fact.
complex<double> TAU[Nw*Nw];
zgeqrf_(&Nw, &Nw, vr_tr, &Nw, TAU, &wkopt, &lwork, &info);
lwork = (int)wkopt.real();
work = new complex<double> [lwork];
zgeqrf_(&Nw, &Nw, vr_tr, &Nw, TAU, work, &lwork, &info);
zungqr_(&Nw, &Nw, &Nw, vr_tr, &Nw, TAU, &wkopt, &lwork, &info);
lwork = (int)wkopt.real();
work = new complex<double> [lwork];
zungqr_(&Nw, &Nw, &Nw, vr_tr, &Nw, TAU, work, &lwork, &info);

//Checking if vr_tr * vr_tr' = I
complex<double> one = 1,zer=0,res[Nw*Nw]; char Tchar = 'T', Nchar = 'N';
zgemm_( &Nchar, &Tchar, &Nw, &Nw, &Nw, &one, vr_tr, &Nw, vr_tr, &Nw, &zer, res, &Nw );
for(i=0;i<Nw;i++){
    for(j=0;j<Nw;j++){
        cout<<res[i*Nw+j]<<"\t";
    }
    cout<<"\n\n";
}

What I get after running this code isn't the identity matrix as it is supposed to be because it should calculate QR factorization from zgeqrf and obtain Q-matrix from zungqr and Q-matrix has to be orthogonal so Q*Q'=I.

What is wrong with this code?

Alireza
  • 410
  • 3
  • 17
  • 2
    Out of curiosity, why are you bothering with lapack when high-level C++ wrappers exist, like [Armadillo](http://arma.sourceforge.net/docs.html#qr)? – The Quantum Physicist Aug 07 '17 at 11:09
  • 1
    @TheQuantumPhysicist Habits and because people I work with use it. Never tried Armadillo though. Can it do all LAPACK routines? – Alireza Aug 07 '17 at 11:17
  • 1
    Check the documentation in the link I provided. When I was still doing research in academia, it had all the routines I needed and even more. For example, LAPACK doesn't have a generic matrix exponentiation function (unless you want to do it yourself by computing the eigenvalues and assuming you have a symmetric matrix), but Armadillo does have that. You can link Armadillo to any LAPACK implementation you wish, and it'll save you so much time. The only issue with Armadillo is that it doesn't support clusters (AFAIK). Otherwise, you should go for it. – The Quantum Physicist Aug 07 '17 at 11:20
  • @TheQuantumPhysicist Awesome! Thanks for the tip! – Alireza Aug 07 '17 at 11:23

2 Answers2

0

I used zunmqr instead of znugqr and it got fixed.

Although I honestly don't know why zungqr didn't work.

Alireza
  • 410
  • 3
  • 17
0

The reason is, the Q-matrix obtained in zungqr has orthonormal columns (with size MxN ) only. So Q^{H}Q=I rather than QQ^{H}=I

Jame
  • 1
  • 1