Низкое потребление ОЗУ 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", ... )

Другие вопросы по тегам