topic Hi Spencer, in Intel® oneAPI Math Kernel Library

MKL Feast

Banerjee__Debasish — Wed, 04 Apr 2018 13:24:20 GMT

Hi,
    I am trying to diagonalize a sparse matrix of dimension NxN, where N = 98310, to find all the eigenvalues and eigenvectors. The matrix is
very sparse (has only 1 in 10^4 non-zero elements), and is symmetric. I store the matrix in sparse CSR format, and use DFEAST_SCSREV
routine (via the Intel MKL library). In the binary form, the matrix only occupies 4MB of space.

    Unfortunately, the program crashes, and the error message "Segmentation fault (core dumped) " is obtained. This output is irrespective of whether I want say 100 eigenvalues in a given interval, or all the eigenvalues. My example code is a minimal modification of the example that intel has.
Could anyone please illuminate the problem? My code is attached below.

with best regards,
Debasish

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mkl_solvers_ee.h"
#include <fstream>
#define max(a, b) (a) < (b) ? (b): (a)

int main()
{
    char          UPLO = 'F'; /* Type of matrix: (F means full matrix, L/U - lower/upper triangular part of matrix) */
    /* Matrix A of size N in CSR format. Size N and 3 arrays are used to store matrix in CSR format */
    
    unsigned int k;
    /* read matrix A from a binary file */
    std::ifstream inFile ("sparse_mat.bin", std::ios::in | std::ios::binary); 
    int nrows, ncols;
    inFile.read((char*)&nrows,sizeof(int));
    inFile.read((char*)&ncols,sizeof(int));
    std::cout<<"No of rows = "<<nrows<<std::endl;
    std::cout<<"No of cols = "<<ncols<<std::endl;

    const MKL_INT N = nrows-1;
    MKL_INT *rows,*cols;
    double *val;
    rows = (MKL_INT*)(calloc(nrows,sizeof(MKL_INT)));
    cols = (MKL_INT*)(calloc(ncols,sizeof(MKL_INT)));
    val =  (double*)(calloc(ncols,sizeof(double)));
    inFile.read((char*)rows, nrows*sizeof(MKL_INT));
    inFile.read((char*)cols, ncols*sizeof(MKL_INT));

    for(k=0;k<ncols;k++) val=-1.0;

    printf("row[10]=%d\n",rows[10]);
    printf("column[10]=%d\n",cols[10]);

    /* Declaration of FEAST variables */
    MKL_INT      fpm[128];      /* Array to pass parameters to Intel MKL Extended Eigensolvers */
    double       Emin, Emax;    /* Lower/upper bound of search interval [Emin,Emax] */

    double       epsout;        /* Relative error on the trace */
    MKL_INT      loop;          /* Number of refinement loop */
    //MKL_INT      L = N;
    MKL_INT      L = 100;
    MKL_INT      M0;            /* Initial guess for subspace dimension to be used */
    MKL_INT      M;             /* Total number of eigenvalues found in the interval */
    double   *E;         /* Eigenvalues */
    double   *X;       /* Eigenvectors */
    double   *res;       /* Residual */
    E = (double*)(calloc(N,sizeof(double)));
    X = (double*)(calloc(N*N,sizeof(double)));
    res = (double*)(calloc(N,sizeof(double)));


    /* Declaration of local variables */
    MKL_INT      info;          /* Errors */
    double       one = 1.0;     /* alpha parameter for GEMM */
    double       zero = 0.0;    /* beta  parameter for GEMM */
    MKL_INT      i, j;
    double       trace, smax, eigabs;

    printf("\n FEAST DFEAST_SCSREV AND DFEAST_SCSRGV EXAMPLE\n");
    /* Initialize matrix X */
    for (i=0; i<N*N; i++){
        X = zero;
    }

    printf("Sparse matrix size %i\n", (int)N);

    /* Search interval [Emin,Emax] */
    Emin = -12.0;
    Emax = 12.0;
    printf("Search interval [ %.15e, %.15e  ]  \n", Emin, Emax);

    M0   = L;
    M    = L;
    loop = 3;
    info = 0;
    epsout = 0.0;

    /* Step 1. Call  FEASTINIT to define the default values for the input FEAST parameters */
    feastinit(
        fpm /* OUT: Array is used to pass parameters to Intel MKL Extended Eigensolvers */
        );

    fpm[0] =  1; /* Extended Eigensolver routines print runtime status to the screen. */

    /* Step 2. Solve the standard Ax = ex eigenvalue problem. */
    printf("Testing dfeast_scsrev routine:\n");
    dfeast_scsrev(
        &UPLO,   /* IN: UPLO = 'F', stores the full matrix */
        &N,      /* IN: Size of the problem */
        val,     /* IN: CSR matrix A, values of non-zero elements */
        rows,    /* IN: CSR matrix A, index of the first non-zero element in row */
        cols,    /* IN: CSR matrix A, columns indices for each non-zero element */
        fpm,     /* IN/OUT: Array is used to pass parameters to Intel MKL Extended Eigensolvers */
        &epsout, /* OUT: Relative error of on the trace */
        &loop,   /* OUT: Contains the number of refinement loop executed */
        &Emin,   /* IN: Lower bound of search interval */
        &Emax,   /* IN: Upper bound of search interval */
        &M0,     /* IN: The initial guess for subspace dimension to be used. */
        E,       /* OUT: The first M entries of Eigenvalues */
        X,       /* IN/OUT: The first M entries of Eigenvectors */
        &M,      /* OUT: The total number of eigenvalues found in the interval */
        res,     /* OUT: The first M components contain the relative residual vector */
        &info    /* OUT: Error code */
        );
    printf("FEAST OUTPUT INFO %d \n",info);
    if ( info != 0 )
    {
        printf("Routine dfeast_scsrev returns code of ERROR: %i", (int)info);
        return 1;
    }
    printf("Number of eigenvalues found %d \n", M);
    //for (i=0; i<M; i++){
    //    printf("%.15e \n", E);
    //}
    
    free(rows); free(cols);
    return 0;
}

Debasish, could you also

Gennady_F_Intel — Wed, 04 Apr 2018 19:04:42 GMT

Debasish, could you also share sparse_mat.bin to investigate this case?

Hi Gennady,

Banerjee__Debasish — Wed, 04 Apr 2018 21:31:00 GMT

Hi Gennady,
I am attaching the file sparse_mat.bin, as sparse_mat.bin_.gz (since in files cannot be uploaded). The bin file has been generated by a code compiled using icpc (intel c++ compiler version 2018) on a Intel(R) Core(TM) i7-4771 CPU @ 3.50GHz processor. Please let me know if something more is needed.
Thank you very much for the help!
Best regards,
Debasish

Hi Debasish, thanks for this

Gennady_F_Intel — Fri, 06 Apr 2018 10:59:54 GMT

Hi Debasish, thanks for this case. We will get back to you soon with any update on this.

Hi,

Banerjee__Debasish — Thu, 21 Jun 2018 16:08:33 GMT

Hi,

Just want to update you all that I don't think Intel FEAST can solve matrices, much beyond N=10000, especially when all
eigenvalues and eigenvectors are needed. Intel MKL lapack routines are much better for this purpose.

best regards,
Debasish

Hi Debasish,

Spencer_P_Intel — Fri, 22 Jun 2018 16:55:00 GMT

Hi Debasish,

A few comments. This matrix you have created has ~98000 egienvalues between -12.75 and 12.75. They are quite close together and so when you are asking for 100 eigenvalues between -12 and 12, if reality you are asking for ~ 98000 eigenvalues. This is not what feast is designed to accomplish. To be more precise, let me explain a little of what feast is doing under the hood. For precise details see https://software.intel.com/en-us/mkl-developer-reference-c-the-feast-algorithm

Feast takes a symmetric sparse matrix and an interval with a decent estimate of how many eigenvalues are to be found on that interval and it creates a (hopefully much smaller) dense matrix that approximates the action of the full matrix on the eigenspace related to the eigenvalues that can be found in that interval. That is, we take our sparse matrix and based on how many eigenvalues are in that interval, we create a dense matrix that can then be used to find those eigenvalues using standard LAPACK dense eigenvalue solvers. The method to do this is outlined in a number of papers and much more detailed than could or should be explained here. It is a rather interesting read if you have the time. (E. Polizzi, Density-Matrix-Based Algorithms for Solving Eigenvalue Problems, Phys. Rev. B. Vol. 79, 115112 (2009) is the main first reference if you are interested)

So, we must pick an interval and get a good estimate of how many eigenvalues are in that interval ( L in the code) We then create a slightly larger subspace size (M0 in the code) and set up our arrays to support M0 eigenvalues being found and returned. I have included here a revised version of your code (called test_feast.cc) that will find the 7 smallest eigenvalues (found between -15 and -10.85 ) The tricky part of using FEAST is getting that interval and estimate. To help in this regard, we have recently added a new functionality to mkl in the 2019 Beta release which allows you to find the N largest or N smallest eigenvalues (smallest here is not smallest magnitude but smallest ordinal) I have also included an example (test_ev.cc) that uses this approach to find smallest 10 eigenvalues of the matrix.

To summarize, you must be very careful in what you ask for when using FEAST. 98k is a large enough number that if you are using 32bit integers, (98k)^2 will actually overflow and so dense LAPACK methods may be able to handle this large if you use the right interfaces (ILP64) and integer types but even if so, it will take a very long time ( O(N^3 ) operations with N~98k is a lot).

Best,

Spencer

Hi Spencer,

Banerjee__Debasish — Fri, 22 Jun 2018 19:48:23 GMT

Hi Spencer,
Thank you so much for your reply, and for providing the improved text example.
It will help me a lot to understand the functionality of the FEAST routine properly.

best,
Debasish