Search code examples
c++matrixlapack

LAPACK QR factorization


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?


Solution

  • I used zunmqr instead of znugqr and it got fixed.

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