【发布时间】：2017-08-08 02:43:11
【问题描述】：
我正在尝试从 LAPACK zgeqrf 和 zungqr 例程中获取 Q-Matrix。
我有一个 Nw×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-matrix 并且 Q-matrix 必须是正交的，所以Q*Q'=I.
这段代码有什么问题？
【问题讨论】：
-
出于好奇，当存在高级 C++ 包装器（如 Armadillo）时，您为什么还要使用 lapack？
-
@TheQuantumPhysicist 习惯和我一起工作的人使用它。虽然从未尝试过犰狳。它可以完成所有 LAPACK 例程吗？
-
检查我提供的链接中的文档。当我还在学术界做研究时，它有我需要的所有例程，甚至更多。例如，LAPACK 没有通用的矩阵求幂函数（除非您想自己通过计算特征值并假设您有一个对称矩阵），但犰狳确实有。您可以将 Armadillo 链接到您希望的任何 LAPACK 实现，它会为您节省大量时间。 Armadillo 的唯一问题是它不支持集群 (AFAIK)。否则，你应该去。
-
@TheQuantumPhysicist 太棒了！感谢您的提示！