Intel MKL LAPACKE_dsyevr на больших матрицах NxN с N~8*10^4

Question

Я пытаюсь диагонализировать (все собственные значения и собственные векторы) матрицы NxN с размерностью N = 80000. Я использую intel MKL 2018 на машине с ~3 ТБ ОЗУ и использую подпрограмму intel LAPACKE_dsyevr. Я прилагаю код ниже, который создает матрицу размера NxN, а затем передает ее в процедуру dsyevr, чтобы найти собственные значения. Хотя код работает для N до ~ 50000, это не работает в случае N ~ 80000, и возникает исключение памяти. Это не может быть связано с памятью моей машины, так как код долго остается в подпрограмме dsyevr до сбоя с ошибкой сегментации.

Обратите внимание, что мне нужно объявить матрицу как std::vector из-за моей основной структуры кода. Поскольку подпрограмма MKL C нуждается в матрице в качестве указателя, я указываю указатель на первый элемент вектора.

Код здесь

    // this uses relatively robust representations 
    #include<iostream>
    #include<fstream>
    #include<stdlib.h>
    #include<stdio.h>
    #include<vector>
    #include<iterator>
    #include<math.h>
    #include "mkl_lapacke.h"
    #include "mkl.h"
    /* Auxiliary routines prototypes */
    extern void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a,   MKL_INT lda );
    extern void fileprint_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda );

    /* Main program */
    int main() {
    unsigned int k;
    unsigned int i,j;
    MKL_INT N, NSQ, LDA, LDZ, NSELECT, info;    
    MKL_INT il, iu, m;
    double abstol, vl, vu;
    N=80000;
    NSQ=N*N;
    LDA=N;
    LDZ=N;

    /* allocate space for the arrays */
    std::vector<MKL_INT> isuppz(2*N,0); 
    std::vector<double> w(N,0.0);
    std::vector<double> z(NSQ,0.0);
    std::vector<double> a(NSQ,0.0);
    std::cout<<"size of vector "<<a.size()<<std::endl;
    std::cout<<"maximum capacity of vector "<<a.max_size()<<std::endl;
    // make the matrix
    k=0;
    for(size_t pos=0; pos < a.size(); pos++){
      i=k/N; j=k%N; 
      a.at(pos) = 1.0/(i+j+1) + (i+j)*exp(-0.344*(i+j));
      k++;
    }
    double* A = &a[0];
    /* print full matrix */
    //print_matrix("Full matrix", N, N, A, LDA);
    /* Executable statements */
    printf( "LAPACKE_dsyevr (column-major, high-level) RRR Example Program  Results\n" );
    /* Solve eigenproblem */
    abstol=-1;
    MKL_INT* ISUPPZ = &isuppz[0];
    double* W = &w[0];
    double* Z = &z[0];
    info = LAPACKE_dsyevr( LAPACK_COL_MAJOR, 'V', 'A', 'U', N, A, LDA,
                        vl, vu, il, iu, abstol, &m, W, Z, LDZ, ISUPPZ );
    printf("Back from the routine and working properly. Info = %d.  Eval[0]=%e\n",info,W[0]);
    //printf("No of eigenvalues = %d\n",m);
    /* Check for convergence */
   if( info > 0 ) {
       printf( "The algorithm failed to compute eigenvalues.\n" );
       exit( 1 );
   }
   isuppz.clear(); 
   w.clear();
   z.clear();
   a.clear();
   exit( 0 );  
   } 

Answer 1

Хорошо, мой друг. Что вы ожидаете? Позвольте мне начать с поздравления с интересной проблемой. Я хотел бы решить это сам. Было сказано, что:

a является 80000*80000*8 = 47.7GB, Это только для хранения вашей матрицы. Но это не все. У вас также есть z и некоторые смехотворно маленькие векторы 80000 двойников. Итак, вы начинаете с выделения 96 ГБ. Сколько оперативной памяти у вас есть для вашего процесса? Но это не твоя проблема, кажется. Ваша проблема LWORK и с этим W. Документация просит LWORK нужно хотя бы 1.26*N:

LWORK is INTEGER
  The dimension of the array WORK.  LWORK >= max(1,26*N).
  For optimal efficiency, LWORK >= (NB+6)*N,
  where NB is the max of the blocksize for DSYTRD and DORMTR
  returned by ILAENV.

  If LWORK = -1, then a workspace query is assumed; the routine
  only calculates the optimal size of the WORK array, returns
  this value as the first entry of the WORK array, and no error
  message related to LWORK is issued by XERBLA.

Но на самом деле, что вы должны сделать, это получить расчет для LWORK из алгоритмов, установив его один раз -1 как предложено выше. Это даст вам правильный номер для LWORK в первой клетке WORK:

Что-то вроде этого:

LWORK=-1;
std::vector<double> w(1);
// Workspace estimation
info = LAPACKE_dsyevr( LAPACK_COL_MAJOR, 'V', 'A', 'U', N, A, LDA,
                    vl, vu, il, iu, abstol, &m, W, Z, LDZ, ISUPPZ );
LWORK = std::static_cast<int>(work.front());
try {
   work.resize(LWORK);
} catch (const std::exception& e) {
   std::cerr << "Failed to allocate work space." << std::endl;
   exit 1;
}
// Actual calculation
info = LAPACKE_dsyevr( LAPACK_COL_MAJOR, 'V', 'A', 'U', N, A, LDA,
                    vl, vu, il, iu, abstol, &m, W, Z, LDZ, ISUPPZ );

Полностью забыл об одной из самых важных вещей, которые я хотел упомянуть: убедитесь, что вам действительно нужна двойная точность. LAPACKE_fsyevr должен иметь оперативную память и быть примерно в два раза быстрее.