topic Found another issue with a in Intel® oneAPI Math Kernel Library

Sparse FEAST

Aaron_M_ — Mon, 09 Sep 2013 17:40:07 GMT

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");
}

}

Aaron,

Zhang_Z_Intel — Mon, 09 Sep 2013 17:58:13 GMT

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?

Zhang,

Aaron_M_ — Mon, 09 Sep 2013 18:03:17 GMT

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.

Aaron,

Zhang_Z_Intel — Tue, 10 Sep 2013 18:05:26 GMT

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 );
    }

}

Zhang,

Aaron_M_ — Wed, 11 Sep 2013 13:28:16 GMT

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.

Zhang,

Aaron_M_ — Wed, 11 Sep 2013 13:33:17 GMT

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.

Found another issue with a

Aaron_M_ — Wed, 11 Sep 2013 15:18:38 GMT

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

]

Aaron,

Zhang_Z_Intel — Wed, 11 Sep 2013 16:28:51 GMT

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!

Here they are. The first one

Aaron_M_ — Wed, 11 Sep 2013 16:34:19 GMT

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

Figured it out! Used the

Aaron_M_ — Wed, 11 Sep 2013 20:49:48 GMT

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

Aaron.

Quote:Aaron M. wrote:

Zhang_Z_Intel — Thu, 12 Sep 2013 16:38:15 GMT

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.

Quick question:

Aaron_M_ — Mon, 16 Sep 2013 20:46:38 GMT

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??

Aaron,

Zhang_Z_Intel — Tue, 17 Sep 2013 16:00:33 GMT

Aaron,

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

Memory is 192GB.

Aaron_M_ — Tue, 17 Sep 2013 17:33:44 GMT

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.

Aaron,

Zhang_Z_Intel — Tue, 17 Sep 2013 17:50:13 GMT

Aaron,

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

I've attached the test file.

Aaron_M_ — Tue, 17 Sep 2013 19:05:54 GMT

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.

Aaron,

Zhang_Z_Intel — Wed, 18 Sep 2013 21:12:44 GMT

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.

Zhang,

Aaron_M_ — Mon, 23 Sep 2013 15:58:03 GMT

Zhang,

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

Aaron.

Aaron,

Zhang_Z_Intel — Mon, 23 Sep 2013 16:39:45 GMT

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).

Aaron,

Gennady_F_Intel — Wed, 30 Oct 2013 06:49:13 GMT

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