Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
This community is designed for sharing of public information. Please do not share Intel or third-party confidential information here.
6587 Discussions

Sparse FEAST

Beginner
592 Views

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

}

19 Replies
Employee
592 Views

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?

Beginner
592 Views

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.

Aaron.

Employee
592 Views

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

}

Beginner
592 Views

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.

Beginner
592 Views

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.

Beginner
592 Views

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

]

Employee
592 Views

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!

Beginner
592 Views

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

Beginner
592 Views

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

Aaron.

Employee
592 Views

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.

Beginner
592 Views

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

Employee
592 Views

Aaron,

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

Beginner
592 Views

Memory is 192GB.

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

Employee
592 Views

Aaron,

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

Beginner
592 Views

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.

Employee
592 Views

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.

Beginner
592 Views

Zhang,

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

Aaron.

Employee
592 Views

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

Moderator
592 Views

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.