5

我是C++ 编程的新手,但我的任务是为非常大的大小矩阵计算对称矩阵(和厄米特矩阵)的特征值和特征向量(标准特征问题Ax=lx ):二项式(L,L/2)其中L大约是18-22。现在我在有大约 7.7 GB 内存的机器上对其进行测试,但最终我可以使用 64 GB 内存的 PC。

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

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

所以我试图将 library 更改为Armadillo library。它工作得很好,但它不如用 all替换的Lapack++好,但当然支持 Hermitian 矩阵。mat Aeigenvec

L=14 的一些统计数据

  • Lapack++ RAM 126MB 时间 7.9s 特征值 + 特征向量

  • 犰狳RAM 216MB 时间 12s 特征值

  • 犰狳RAM 396MB 时间 15s 特征值+特征向量

让我们做一些计算:double变量大约是8B。我的矩阵有大小 binomial(14,7) = 3432,所以在理想情况下它应该有3432^2*8/1024^2 = 89 MB

我的问题是:是否可以修改或强制犰狳做一个像Lapack++一样的好把戏?犰狳的用途LAPACKBLAS惯例。或者也许有人可以推荐使用另一个库来解决这个问题的另一种方法?

PS:我的矩阵真的很稀疏。它有大约2 * binomial(L,L/2)个非零元素。我尝试用SuperLUCSC 格式计算它,但它不是很有效,因为 L=14 -> RAM 185MB,但时间为 135 秒。

4

2 回答 2

7

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)。如果你想保持密集的矩阵格式,你也可以试试ComplexEigenSolverEigen 库的。

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 示例。

于 2015-08-28T22:52:49.430 回答
1

如果将来有人遇到与我相同的问题,在我找到解决方案后会有一些提示(感谢所有发布一些答案或线索的人!)。

在英特尔网站上,您可以找到一些用 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( ... )to zheev_( ... )(只需通过替换就可以在代码中使用它,但我不知道它为什么起作用。我通过一些实验弄清楚了。)

或者,您可以将printf函数转换为标准流std::cout并将包含的标头替换stdio.hiostream.

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

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

于 2015-08-30T09:27:00.817 回答