броненосец — низкий потребляющий ОЗУ c ++ собственный решатель

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

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

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

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

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

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

Реальный вопрос заключается в том, нужны ли вам все собственные векторы, все собственные значения или только самые большие собственные значения. Потому что в последнем случае есть действительно эффективные методы. Посмотрите на Арнольди /Lanczos Итеративные алгоритмы. Огромный прирост памяти возможен, если матрица разрежена, поскольку выполняются только матрично-векторные произведения: нет необходимости сохранять плотный формат. Это то, что делается в библиотеке SlepC, которая использует форматы разреженной матрицы Petsc. Вот пример Slepc который может быть использован в качестве отправной точки.