Sparse FEAST

[Aaron_M_]

Hi,

I ran into an issue with the sparse FEAST implementation of "dfeast_scsrev" & "dfeast_scsrgv".
The actual matrix I need eigen-decomposed is of size 4096 x 4096. Here is a test case that
presents the same problem:

A = [ 1 1 0
      1 0 0
      0 0 0]

In CSR format A is:
vals=(1, 1, 1)
rows=(1, 3, 4)
cols=(1, 2, 1).

Now this does not take in to accoun the third row and column of all zeros. In fact, in CSR format
it is no different than:

A'= [ 1 1
      1 0].

However, for my problem, I can't simply truncate the Hilbert space, it is imperative that I return
an eigenvalue of 0. The dense version of FEAST( dfeast_syev) can eigen-decompose both these matrices,
however, the sparse version can only do the latter, not the former. It just hangs. Please let me know
of any suggestions! Thank you!

Here is the source code for the general version of Sparse:
Version: icc version 13.1.0 (gcc version 4.4.7 compatibility)
///////////Trial for Sparse General///////////

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include <vector>
#include <vector>

using namespace std;

void convertDense2CSR( double *A, int N, vector<double> &Values, vector<int> &Rows, vector<int> &Cols) {
    // Convert matrix A (dense format) to CSR

    Values.clear();
    Rows.clear();
    Cols.clear();

    int nnz = 0;
    int lastRow = -1;
    for (int i=0; i<N; i++) {
        for (int j=0; j<N; j++) {
            if (A[i*N+j] != 0) {
                Values.push_back( A[i*N+j] );
                Cols.push_back( j+1 );
                if (i != lastRow) {
                    Rows.push_back( nnz+1 );
                }
                nnz++;
                lastRow = i;
            }
        }
    }
    Rows.push_back( nnz+1 );
}

int main()
{

  int N = 3;
  int i, j;
  double A[9] = {1, 1, 0, 1, 0, 0, 0, 0, 0};   
  double B[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
    printf (" 'A: \n\n");
    for (i=0; i<3; i++) {
    for (j=0; j<3; j++) {
    printf ("%12.4f", A[i*N+j]);
    }
    printf ("\n");
    }
    printf (" B: \n\n");
    for (i=0; i<3; i++) {
    for (j=0; j<3; j++) {
    printf ("%12.4f", B[i*N+j]);
    }
    printf ("\n");
    }

  ////////////////////////////////
  /////////////FEAST//////////////
  ////////////////////////////////

   
  /* EIGEN-VALUE SYSTEM*/

  int LDA = N;
  char UPLO = 'F';
  double Emin = -10.0, Emax = 10.0;
  int M0 = N;

  /* INPUT PARAMETERS FOR FEAST*/
  int feastparam[128];
 
  /* OUTPUT VARIABLES FOR FEAST*/
  double  *E, *res, *X;
  double  epsout, trace;
  int  loop, info, M;

  /*Allocate memory for eigenvaluues. eigenvectors/residuals*/
  /* Note: I'm doing this differently than the example*/
      
 
  E = (double *)calloc( M0*sizeof( double ), 64 );  //eigenvalues
  res = (double *)calloc( M0*sizeof( double ), 64 );  //residuals
  X = (double *)calloc( M0*M0*sizeof( double ), 64 ); //eigenvectors

  /* !!!!!!!!!!!!! FEAST !!!!!!!!!!!!!!!!! */
  FEASTINIT(feastparam);
  feastparam[0] = 0; // change default value
  //feastparam[5] = 0; // change default value
  dfeast_syev( &UPLO, &N, A, &LDA, feastparam, &epsout, &loop, &Emin, &Emax, &M0, E, X, &M, res, &info );
  ///*
  printf("\nReport using Dense version:\n");
    printf(" Eigenvalues\n\n");
    for (i=0; i<N; i++) {
    printf ("%12.5f \n", E);
    }
 
    printf("\n");//*/
    ///*
 
 
    printf (" 'X' Eigenvectors \n\n");
    for (i=0; i<N; i++) {
    for (j=0; j<N; j++) {
    printf ("%12.4f", X[i*N+j]);
    }
    printf ("\n");
    }
    //*/
    printf("\nReport using General  Sparse version:\n");

    //char UPLO = 'F';
  //double Emin = -10.0, Emax = 10.0;
  int M01 = N;

  /* INPUT PARAMETERS FOR FEAST*/
  //int feastparam[128];
 
  /* OUTPUT VARIABLES FOR FEAST*/
  double  *E1, *res1, *X1;
  double  epsout1, trace1;
  int  loop1, info1, M1;

  /*Allocate memory for eigenvaluues. eigenvectors/residuals*/


  vector<double> Values;
  vector<int>  Cols, Rows;
  vector<double> ValuesB;
  vector<int>  ColsB, RowsB;

  convertDense2CSR( A, N, Values, Rows, Cols );
  convertDense2CSR( B, N, ValuesB, RowsB, ColsB );

 
  /*
  cout << "Set CSR version of A.  We have, values: " << endl;
  for (int i=0; i< Values.size(); i++) {
    cout << Values.at(i) << "    ";
  } cout << endl;
  cout << "We have, cols: " << endl;
  for (int i=0; i< Cols.size(); i++) {
    cout << Cols.at(i) << "    ";
  } cout << endl;
  cout << "We have, rows: " << endl;
  for (int i=0; i< Rows.size(); i++) {
    cout << Rows.at(i) << "    ";
  } cout << endl;
  */
 

  E1 = (double *)calloc( M01*sizeof( double ), 64 );  //eigenvalues
  res1 = (double *)calloc( M01*sizeof( double ), 64 );  //residuals
  X1 = (double *)calloc( M01*M01*sizeof( double ), 64 ); //eigenvectors

  /* !!!!!!!!!!!!! FEAST !!!!!!!!!!!!!!!!! */
  FEASTINIT(feastparam);
  feastparam[0] = 0; // change default value
 
 

  printf("FEAST OUTPUT INFO %d \n",info);
  dfeast_scsrgv(&UPLO,&N,(&Values[0]),(&Rows[0]),(&Cols[0]),(&ValuesB[0]),(&RowsB[0]),(&ColsB[0]) ,feastparam,&epsout1,&loop1,&Emin,&Emax,&M01,E1,X1,&M1,res1,&info1);
  if ( info1 != 0 ){ return 1;

  }

  printf(" EigenValues2: \n");

  for (i=0; i<N; i++) {
    printf ("%12.5f", E1);
  }
    
  printf("\n Eigenvector matrix X2: \n");

  for (i=0; i<N; i++) {
    for (j=0; j<N; j++) {
      printf ("%12.4f", X1[j+i*N]);
    }
    printf("\n");
  }

}

[Zhang_Z_Intel]

Aaron,

Thanks for your question. I'll do some investigation and let you know. For the time being, can you use dense FEAST (dfeast_syev) as workaround?

[Aaron_M_]

Zhang,

As of right now I am. But it is a huge bottleneck in my program, hence the allure of using the sparse version. My 4096x4096 is very sparse, it only has 600 non-zero elements in it, but it takes a few mins for the dense FEAST run. I was hoping that by using sparse I could shave that down to under a minute.

Thanks for the reply,

Aaron.

[Zhang_Z_Intel]

Aaron,

It turns out that there is an error in your CSR representation of A. CSR format stipulates that the length of 'rows' must be (number_of_rows + 1). So 'rows' should be {1, 3, 4, 4} in your code. I've changed the convertDense2CSR routine accordingly and now your test code works fine. See below:

void convertDense2CSR( double *A, int N, vector<double> &Values, vector<int> &Rows, vector<int> &Cols) {
    // Convert matrix A (dense format) to CSR

    Values.clear();
    Rows.clear();
    Cols.clear();

    int nnz = 0;
    int lastRow = -1;
    for (int i=0; i<N; i++) {
        for (int j=0; j<N; j++) {
            if (A[i*N+j] != 0) {
                Values.push_back( A[i*N+j] );
                Cols.push_back( j+1 );
                if (i != lastRow) {
                    Rows.push_back( nnz+1 );
                }
                nnz++;
                lastRow = i;
            }
        }
    }
    Rows.push_back( nnz+1 );
  
    if (Rows.size() < N+1) {
        Rows.push_back( nnz+1 );
    }
  
}

[Aaron_M_]

Zhang,

Thanks so much for the help. I actually have tried to implement this into another code, with a matrix of size 64x64. Now instead of hanging, I get an info output of -4. There isn't any documentation on that and it isn't listed as one of the errors. Any suggestions?

Aaron.

[Aaron_M_]

Zhang,

Thanks so much for the help! I've tried implementing the change in a larger state space, 64x64 matrix. The FEAST info output is -4. That is not listed in the documentation and was wondering if you had seen something like that before?

Aaron.

[Aaron_M_]

Found another issue with a test case. Now, unfortunately this still doesn't resolve my issue but thought someone out there may have seen this:

if A=[3 0 0 0; 0 0 0 0 ; 0 0 0 0; 0 0 0 5]

then although the sparse feast output is 0, there is no output from the eigen-vectors or values. I've attached my code if someone wouldn't mind taking a look. I put a "while loop" in but I get the wrong eigen values and vectors although I get an output. The while-loop is in place of the if statement in my function to convert Dense2Csr. It is commented out.

000

]

[Zhang_Z_Intel]

Aaron,

Please attach your test code for the case A=[3 0 0 0; 0 0 0 0 ; 0 0 0 0; 0 0 0 5]. Thanks!

[Aaron_M_]

Here they are. The first one "sparse2" is the 4 x 4. The second one is a 64 x 64.

[Aaron_M_]

Figured it out! Used the mklddnscsr function instead of use my own. It works like a charm! Thanks guys!

Aaron.

[Zhang_Z_Intel]

Aaron M. wrote:

Figured it out! Used the mklddnscsr function instead of use my own. It works like a charm! Thanks guys!

Aaron.

Fantastic! Thank you for letting us know.

[Aaron_M_]

Quick question:

I have two implementations of FEAST(dense and sparse) of a 4096x4096 matrix. The dense version works great, the sparse version has an output that I attached. It can't be because the data set is too large because the dense version works, and the mklddncsr works great too. Any ideas??

[Zhang_Z_Intel]

Aaron,

How many OpenMP threads did you specify when you ran the sparse FEAST? How big is the memory on your system?

[Aaron_M_]

Memory is 192GB.

OMP_NUM_THREADS is 16.

Initially I thought that there just wasn't enough memory, but it doesn't seem to be a problem with dense FEAST.

[Zhang_Z_Intel]

Aaron,

This looks like a bug in sparse FEAST. Can you share with us your test matrix? Thanks!

[Aaron_M_]

I've attached the test file. Unfortunately I have to jump through hoops to get to the test matrix I need eigen-decomposed. I'm sorry, it's un-inevitable. The FEAST implementation starts at line 595 and is marked "FEAST". From that section on:

1. I use mklddnscr to convert my test matrix labeled "Sym" to "Symcsr , Sym J, & Sym I"

2.  at line 662 is where I call dfeast_scsrev. It will crash. 

3. line 669 is a commented out dfeast_syev. That works. Just comment out line 662 and un-comment dense FEAST and the code runs beautifully, albiet slowly. It takes ~140 seconds to work with dense FEAST.

Thank you.

[Zhang_Z_Intel]

Aaron,

The MKL engineering team confirms this is a bug. Thank you for reporting it. I will keep you updated when the big is fixed.

[Aaron_M_]

Zhang,

Has there been any progress? Please let me know. Thank you.

Aaron.

[Zhang_Z_Intel]

Aaron,

The MKL team is looking at this issue. You will hear back from us within 2 weeks. We will try to provide a temporary workaround, if possible. A permanent fix to this problem will be put in the product in the next update release (due in 2~3 months).

[Gennady_F_Intel]

Aaron, 

please check if the problem is still exists in the latest update 1 ( MKL v.11.1 Update 1) released the last Friday and let us know the results.

Gennady