Низкое потребление ОЗУ C++
Я новичок в программировании на C++, но у меня есть задача вычислить собственные значения и собственные векторы (стандартная проблема Ax = lx) для симметричных матриц (и эрмитовых)) для очень большой матрицы размера: биномиальная (L,L/2), где L около 18-22. Сейчас я тестирую его на машине с 7,7 ГБ оперативной памяти, но в итоге у меня будет доступ к ПК с 64 ГБ ОЗУ.
Я начал с Lapack++. В начале мой проект предполагал решить эту проблему только для симметричных вещественных матриц.
Эта библиотека была великолепна. Очень быстрый и маленький объем оперативной памяти. Он имеет возможность вычислить собственные векторы и поместить во входную матрицу A для экономии памяти. Оно работает! Я думал, что Lapack++ eigensolver может работать с эрмитовой матрицей, но не может по неизвестной причине (возможно, я делаю что-то не так). Мой проект развился, и я должен быть в состоянии рассчитать эту проблему для эрмитовых матриц.
Поэтому я попытался изменить библиотеку на библиотеку Armadillo. Он отлично работает, но он не так хорош, как Lapack++, который заменяет mat A
со всем eigenvec
, но, конечно, поддерживает эрмитовы матрицы.
Немного статистики для L = 14
Lapack++ RAM 126MB время 7,9 с собственное значение + собственные векторы
Armadillo RAM 216MB время 12с собственное
Armadillo RAM 396MB время 15с собственное значение + собственные векторы
Давайте сделаем некоторые расчеты: double
переменная составляет около 8B. Моя матрица имеет размер binomial (14,7) = 3432, поэтому в идеальном случае она должна иметь 3432^2*8/1024^2 = 89 МБ.
Мой вопрос: возможно ли изменить или проинформировать Армадилло, чтобы он сделал хороший трюк, как Lapack++? Армадилло использует LAPACK
а также BLAS
Подпрограммы. Или, может быть, кто-то может порекомендовать другой подход к этой проблеме, используя другую библиотеку?
PS: моя матрица действительно разреженная. Он имеет около 2 * биномиальных (L,L/2) ненулевых элементов. Я пытался вычислить это с SuperLU
в формате CSC, но это было не очень эффективно, для L=14 -> RAM 185MB, но время 135 с.
2 ответа
И Лапакпп, и Армадилло полагаются на Лапака для вычисления собственных значений и собственных векторов комплексных матриц. Библиотека Лапака предоставляет различные способы выполнения этих операций для сложных эрмитовых матриц.
Функция
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 также могут быть эффективны в этом аспекте.
Реальный вопрос заключается в том, нужны ли вам все собственные векторы, все собственные значения или только самые большие собственные значения. Потому что в последнем случае есть действительно эффективные методы. Взгляните на итерационные алгоритмы Арнольди / Ланцоша. Огромный прирост памяти возможен, если матрица является разреженной, поскольку выполняются только матрично-векторные произведения: нет необходимости сохранять плотный формат. Это то, что делается в библиотеке SlepC, которая использует форматы разреженной матрицы Petsc. Вот пример Slepc, который можно использовать в качестве отправной точки.
Если у кого-то возникнет та же проблема, что и у меня в будущем, то после того, как я нашел решение, есть несколько советов (спасибо за все, что опубликовали некоторые ответы или подсказки!).
На сайте Intel вы можете найти несколько хороших примеров, написанных на фортране и 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( ... )}
Изменить функцию вызова вызывающего Лапака _
символ например: zheev( ... )
в zheev_( ... )
(сделать это для всех в коде просто путем подстановки, но я не понимаю, почему это работает. Я понял это, проведя некоторые эксперименты.)
При желании вы можете конвертировать printf
функция в стандартный поток std::cout
и заменить включенные заголовки stdio.h
в iostream
,
Чтобы скомпилировать команду run, выполните следующие действия: g++ test.cpp -o test -llapack -lgfortran -lm -Wno-write-strings
О последнем варианте -Wno-write-strings
Теперь я не знаю, что он делает, но, вероятно, есть проблемы с примером Intel, когда они помещают строку, а не символ в вызывающую функцию. zheev( "Vectors", "Lower", ... )