低内存消耗c++特征解算器

Lapackpp和Armadillo都依赖于Lapack来计算复矩阵的特征值和特征向量。Lapack库提供了对复杂厄米矩阵执行这些操作的不同方法。

• 函数zgeev()不关心矩阵是否是厄米矩阵。这个函数由Lapackpp库在函数LaEigSolve中为类型为LaGenMatComplex的矩阵调用。Armadillo库的eig_gen()函数调用这个函数

• 函数zheev()专用于复厄米矩阵。首先调用ZHETRD将厄米矩阵化简为三对角线形式。根据是否需要特征向量，然后使用QR算法计算特征值和特征向量(如果需要)。如果选择了方法std, Armadillo库的eig_sym()函数将调用此函数

• 如果不需要特征向量，函数zheevd()与zheev()做同样的事情。否则，它会使用分治算法(参见zstedc())。如果选择了方法dc, Armadillo库的eig_sym()函数将调用该函数。由于分治法被证明对于大矩阵更快，所以它现在是默认方法。

有更多选项的函数在Lapack库中可用。(参见zheevr()或zheevx)。如果你想保持密集的矩阵格式，你也可以尝试Eigen库的ComplexEigenSolver。

下面是一个使用Lapack库的C包装器LAPACKE的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迭代算法。如果矩阵是空间的，则可能获得巨大的内存增益，因为只执行矩阵-向量乘积:没有必要保持密集格式。这是在SlepC库中完成的，它利用了Petsc的空间矩阵格式。下面是一个可以作为起点的Slepc示例:

如果有人有同样的问题，我在未来有一些提示后，我找到了解决方案(感谢所有张贴的一些答案或线索!)。

在英特尔网站上，你可以找到一些用Fortran和c写的很好的例子。例如厄米特征值问题例程zheev():https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/zheev_ex.c.htm

要使它在c++中工作，你应该在代码中编辑一些行:

在原型函数声明中，对所有函数做类似的声明:extern void zheev( ... )更改为extern "C" {void zheev( ... )}

更改调用lapack的函数，添加_符号，例如:zheev( ... )到zheev_( ... )(通过替换使其在代码中全部可用，但我不知道为什么它可以工作。我通过做一些实验弄明白了。

您可以选择将printf函数转换为标准流std::cout，并将包含的标头stdio.h替换为iostream。

编译运行命令:g++ test.cpp -o test -llapack -lgfortran -lm -Wno-write-strings

关于最后一个选项-Wno-write-strings，我现在不知道它在做什么，但可能有问题英特尔的例子，当他们把字符串而不是字符调用函数zheev( "Vectors", "Lower", ... )