Lapackpp 和 Armadillo 都依赖 Lapack 来计算复杂矩阵的特征值和特征向量。Lapack 库为复杂的厄米矩阵提供了不同的方法来执行这些操作。
该函数zgeev()
不关心矩阵是 Hermitian。该函数由 Lapackpp 库调用,用于LaGenMatComplex
函数中的类型矩阵LaEigSolve
。Armadillo 库的函数eig_gen()
调用此函数。
该函数zheev()
专用于复杂的 Hermitian 矩阵。它首先调用 ZHETRD 将 Hermitian 矩阵简化为三对角形式。根据是否需要特征向量,然后使用QR 算法来计算特征值和特征向量(如果需要)。如果方法被选中,犰狳库的函数eig_sym()
调用这个函数。std
该函数的作用zheevd()
与不需要特征向量相同zheev()
。否则,它会使用分而治之的算法(请参阅 参考资料zstedc()
)。如果方法被选中,犰狳库的函数eig_sym()
调用这个函数。dc
由于分而治之被证明对于大型矩阵更快,因此它现在是默认方法。
Lapack 库中提供了具有更多选项的函数。(见zheevr()
或zheevx
)。如果你想保持密集的矩阵格式,你也可以试试ComplexEigenSolver
Eigen 库的。
LAPACKE
这是一个使用Lapack 库的 C 包装器的小 C++ 测试。它是由g++ main.cpp -o main2 -L /home/...../lapack-3.5.0 -llapacke -llapack -lblas
#include <iostream>
#include <complex>
#include <ctime>
#include <cstring>
#include "lapacke.h"
#undef complex
using namespace std;
int main()
{
//int n = 3432;
int n = 600;
std::complex<double> *matrix=new std::complex<double>[n*n];
memset(matrix, 0, n*n*sizeof(std::complex<double>));
std::complex<double> *matrix2=new std::complex<double>[n*n];
memset(matrix2, 0, n*n*sizeof(std::complex<double>));
std::complex<double> *matrix3=new std::complex<double>[n*n];
memset(matrix3, 0, n*n*sizeof(std::complex<double>));
std::complex<double> *matrix4=new std::complex<double>[n*n];
memset(matrix4, 0, n*n*sizeof(std::complex<double>));
for(int i=0;i<n;i++){
matrix[i*n+i]=42;
matrix2[i*n+i]=42;
matrix3[i*n+i]=42;
matrix4[i*n+i]=42;
}
for(int i=0;i<n-1;i++){
matrix[i*n+(i+1)]=20;
matrix2[i*n+(i+1)]=20;
matrix3[i*n+(i+1)]=20;
matrix4[i*n+(i+1)]=20;
matrix[(i+1)*n+i]=20;
matrix2[(i+1)*n+i]=20;
matrix3[(i+1)*n+i]=20;
matrix4[(i+1)*n+i]=20;
}
double* w=new double[n];//eigenvalues
//the lapack function zheev
clock_t t;
t = clock();
LAPACKE_zheev(LAPACK_COL_MAJOR,'V','U', n,reinterpret_cast< __complex__ double*>(matrix), n, w);
t = clock() - t;
cout<<"zheev : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
cout<<"largest eigenvalue="<<w[n-1]<<endl;
std::complex<double> *wc=new std::complex<double>[n];
std::complex<double> *vl=new std::complex<double>[n*n];
std::complex<double> *vr=new std::complex<double>[n*n];
t = clock();
LAPACKE_zgeev(LAPACK_COL_MAJOR,'V','V', n,reinterpret_cast< __complex__ double*>(matrix2), n, reinterpret_cast< __complex__ double*>(wc),reinterpret_cast< __complex__ double*>(vl),n,reinterpret_cast< __complex__ double*>(vr),n);
t = clock() - t;
cout<<"zgeev : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
cout<<"largest eigenvalue="<<wc[0]<<endl;
t = clock();
LAPACKE_zheevd(LAPACK_COL_MAJOR,'V','U', n,reinterpret_cast< __complex__ double*>(matrix3), n, w);
t = clock() - t;
cout<<"zheevd : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
cout<<"largest eigenvalue="<<w[n-1]<<endl;
t = clock();
LAPACKE_zheevd(LAPACK_COL_MAJOR,'N','U', n,reinterpret_cast< __complex__ double*>(matrix4), n, w);
t = clock() - t;
cout<<"zheevd (no vector) : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
cout<<"largest eigenvalue="<<w[n-1]<<endl;
delete[] w;
delete[] wc;
delete[] vl;
delete[] vr;
delete[] matrix;
delete[] matrix2;
return 0;
}
我的电脑输出是:
zheev : 2.79 seconds
largest eigenvalue=81.9995
zgeev : 10.74 seconds
largest eigenvalue=(77.8421,0)
zheevd : 0.44 seconds
largest eigenvalue=81.9995
zheevd (no vector) : 0.02 seconds
largest eigenvalue=81.9995
这些测试可以通过使用 Armadillo 库来执行。直接调用 Lapack 库可能会让您获得一些内存,但 Lapack 的包装器在这方面也很有效。
真正的问题是您是否需要所有特征向量、所有特征值或仅需要最大特征值。因为在最后一种情况下确实有有效的方法。看看 Arnoldi/ Lanczos迭代算法。如果矩阵是稀疏的,则可能会获得巨大的内存增益,因为只执行矩阵向量乘积:不需要保持密集格式。这就是在使用 Petsc 的稀疏矩阵格式的 SlepC 库中所做的。这是一个可以用作起点的 Slepc 示例。