LAPACK QR因子分解

我试图从LAPACK zgeqrf和zungqr例程中获得Q-Matrix。LAPACK QR因子分解

我有一个Nw-by-Nw复数矩阵，其列上有非正交向量。这是我的C++代码。（矩阵名为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"; 
}

什么运行此代码后，我得到的是不是单位矩阵，因为它应该是，因为它应该从zgeqrf计算QR分解，并从zungqr和Q获得Q-矩阵-matrix必须是正交的，所以Q*Q'=I。

这段代码有什么问题？

来源

2017-08-07 Alireza

出于好奇，当高级C++包装器存在时，为什么你会用lapack打扰，像[Armadillo]（http://arma.sourceforge.net/docs.html#qr）？ –

@TheQuantumPhysicist习惯，因为我工作的人使用它。尽管没尝试过犰狳。它可以执行所有LAPACK例程吗？ – Alireza

检查我提供的链接中的文档。当我还在学术界进行研究时，它有我需要的所有例程，甚至更多。例如，LAPACK没有一个通用的矩阵求幂函数（除非你想通过计算特征值并假设你有一个对称矩阵来完成它），但是Armadillo确实有这个功能。您可以将Armadillo链接到您希望的任何LAPACK实现，并且可以为您节省很多时间。 Armadillo唯一的问题是它不支持集群（AFAIK）。否则，你应该去做。 –

我使用zunmqr而不是znugqr它得到了修复。

虽然我真的不知道为什么zungqr没有工作。

来源

2017-08-07 18:43:11 Alireza

LAPACK QR因子分解

回答

相关问题