собственный решатель С++ с низким потреблением оперативной памяти

И Lapackpp, и Armadillo полагаются на Lapack для вычисления собственных значений и собственных векторов комплексных матриц. Библиотека Лапака предоставляет другой способ выполнения этих операций для сложных эрмитовых матриц.

• Функция zgeev() не заботится о том, чтобы матрица эрмитов. Эта функция вызывается библиотекой Lapackpp для матриц типа LaGenMatComplex в функции LaEigSolve. Эта функция вызывается функцией eig_gen() библиотеки Armadillo.

• Функция zheev() предназначена для сложных эрмитовых матриц. Сначала он вызывает ЖЕТРД, чтобы привести эрмитову матрицу к трехдиагональной форме. В зависимости от того, нужны ли собственные векторы, он использует QR-алгоритм для вычисления собственных значений и собственных векторов. (если нужно). Функция eig_sym() библиотеки Armadillo вызывает эту функцию, если выбран метод std.

• Функция zheevd() делает то же самое, что и zheev(), если собственные векторы не требуются. В противном случае используется алгоритм «разделяй и властвуй» (см. zstedc()). Функция eig_sym() библиотеки Armadillo вызывает эту функцию, если выбран метод dc. Поскольку метод «разделяй и властвуй» оказался быстрее для больших матриц, теперь он является методом по умолчанию.

Функции с большим количеством опций доступны в библиотеке Lapack. (см. zheevr() или zheevx). Если вы хотите сохранить формат плотной матрицы, вы также можете попробовать ComplexEigenSolver библиотеки Eigen.

Вот небольшой тест C++ с использованием C-обертки LAPACKE библиотеки Lapack. Составил 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 который можно использовать в качестве отправной точки.

Если у кого-то в будущем возникнет та же проблема, что и у меня, есть несколько советов после того, как я нашел решение (спасибо всем, кто опубликовал ответы или подсказки!).

На сайте Intel вы можете найти несколько хороших примеров, написанных на Fortran и C. Например, процедура эрмитовой задачи на собственные значения zheev(): https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/zheev_ex.c.htm