低内存消耗的c++特征求解器

Question

我是C++编程中的新手，但我有一项任务要计算对称矩阵(和hermitian)的特征值和特征向量(标准特征问题Ax=lx)，它们的大小很大:二项式(L，L/2)，其中L约为18-22。现在我正在机器上测试它，它有大约7.7GB的内存，但最终我可以访问64 GB内存的PC。

我从Lapack++开始。在开始的时候，我的项目假设只对对称实矩阵解决这个问题。

这个图书馆很棒。非常快和小内存消耗。它可以选择计算特征向量并将其放入输入矩阵A中以节省内存。它起作用了!我认为Lapack++ eigensolver可以处理Hermitian矩阵，但由于未知的原因它不能处理(也许我做错了什么)。我的项目已经发展，我也应该能够计算这个问题的厄米矩阵。

所以我试着把图书馆改成Armadillo图书馆。它工作得很好，但是它并不像Lapack++那样好，它用所有的eigenvec替换了mat A，但是当然支持hermitian矩阵。

关于L=14的一些统计

• Lapack++ RAM 126 RAM时间7.9s特征向量+特征向量
• 鲤鱼RAM 216 RAM时间12s特征值
• 鲤鱼RAM 396 15s时间15s eigenval+eigenvectors

让我们计算一下：double变量大约是8B。我的矩阵的大小为二项式(14，7)= 3432，因此在理想情况下，它应该有3432^2*8/1024^2 = 89 MB。

我的问题是:修改或强制Armadillo作为Lapack++做一个很好的技巧是可行的吗？Armadillo使用LAPACK和BLAS例程。或者有人可以用另一个库来解决这个问题呢？

我的矩阵很稀疏。它有约2*二项式(L，L/2)非零元素。我尝试过用SuperLU以CSC格式计算这个值，但是对于L=14 -> RAM 185 in，但是对于time 135 s，它并不是很有效。

Answer 1

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

• 函数zgeev()并不关心矩阵是Hermitian。函数LaGenMatComplex中的LaEigSolve类型矩阵的Lapackpp库调用此函数。Armadillo库的函数eig_gen()调用此函数。
• 函数zheev()专用于复Hermitian矩阵。它首先调用ZHETRD将Hermitian矩阵简化为三对角形式。根据是否需要特征向量，然后使用QR算法计算特征值和特征向量(如果需要的话)。如果选择了方法eig_sym()，则Armadillo库的函数std调用此函数。
• 如果不需要特征向量，则函数zheevd()执行与zheev()相同的操作。否则，它将使用除法和征服器算法(请参阅zstedc())。如果选择了方法eig_sym()，则Armadillo库的函数dc调用此函数。由于分而治之对于大型矩阵被证明是更快的，所以现在它是默认的方法。

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

下面是使用Lapack库的C包装器LAPACKE进行的一个小小的LAPACKE测试。它是由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/兰索斯独立算法。如果矩阵是稀疏的，则可以获得巨大的内存增益，因为只执行矩阵向量积:不需要保持密集的格式。这就是在SlepC库中所做的工作，它使用了Petsc的稀疏矩阵格式。下面是一个搜索器的例子可以作为一个起点。

Answer 2

如果将来有人和我有同样的问题，在我找到解决方案后有一些建议(谢谢你发布了一些答案或线索！)

在英特尔网站上，您可以找到一些用Fortran和C编写的很好的例子，例如hermitian特征值问题例程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", ... )中添加字符串而不是字符时，英特尔的例子可能有问题