Trouble using MKL DSS subroutines

[vijay1]

I am having trouble using MKL DSS subroutines to solve a symmetric sparse matrix in C++. I called the subroutines that intialize and decompose the matrix in one function say Decompose() function: dss_create(), dss_define_structure(), dss_reorder(), dss_factor_real(). Then I called the solve subroutine: dss_solve_real() in another function say Solve(). The handle is defined as a class variable and is initialized in the Decompose() function. I first execute the Decompose() function and then call Solve(). When the code reaches the subroutine dss_solve_real() it is crashing out. But if I put all the subroutines in one function it works ok. Am I not allowed to call these subroutines from different function? Is there a local memory that is destroyed when I leave the first function? I would greatly appreciate if somebody has any suggestions to fix this problem.
Thanks,
Vijay

[Konstantin_A_Intel]

Hi Vijay,

Could you please send us a small reproducer? You can also post some pieces of the code right here in the forum. I mostly interested in declarations of the functions, of the class which stores handle and calling program.

As it works fine when all the calls made in one place - it seems that something wrong happened with DSS parameters between calling of 2 subroutines.

Regards,
Konstantin

[vijay1]

Hi Konstantin,

Thanks for your reply. I figured out the reason for the crash. Basically I had an existing class with specific data structure to store matrix contents and for testing the MKL DSS I created a new pointer to double with values of the matrix and I was mistakenly deleting this pointer once I exist the function Decompose() (thinking that the dss handle now stores this data). I guess the dss handle stores only the pointer to these values. Do you think this makes sense?

Regards,
Vijay

[vijay1]

Hi Konstantin,

I am now having trouble solving a complex system of equations using MKL DSS. The solution I obtain from these subroutines is completely wrong. Could you please advise as to what I could be doing wrong here? Below is my code for your reference:

Thanks,
Vijay

#include
#include
#include
#include
#include
using namespace std;
#include "mkl_dss.h"
#include "mkl_dss.h"
#include "mkl_types.h"

/*
** Define the array and rhs vectors
*/
#define NROWS 5
#define NCOLS 5
#define NNONZEROS 9
#define NRHS 1
int ROWINDEX[NROWS+1] = { 1, 6, 7, 8, 9, 10 };
int COLUMNS[NNONZEROS] = { 1, 2, 3, 4, 5, 2, 3, 4, 5 };
//double VALUES[NNONZEROS] = { 9, 1.5, 6, .75, 3, 0.5, 12, .625, 16 };
complex VALUES[NNONZEROS] = {complex (9,3), complex (1.5,1), complex (6,2), complex (0.75,1), complex (3,1.2), complex (0.5,0.5), complex (12,9), complex (0.625,1), complex (16,7)};
//double RHS[NCOLS] = { 1, 2, 3, 4, 5 };
complex RHS[NCOLS] = {complex (1,1), complex (2,2), complex (3,1.5), complex (4,1), complex (5,1.2)};

class SymmSquareMatrix{
public:
int nRows;
int nCols;
int nNonZeros;
//complex * values;
_DOUBLE_COMPLEX_t *values;
int* rowIndex;
int* columns;
_MKL_DSS_HANDLE_t handle;
bool Decompose();
bool Solve();
SymmSquareMatrix();
virtual ~SymmSquareMatrix();
};

SymmSquareMatrix::SymmSquareMatrix(){
//define the matrix here!
nRows=NROWS;
nCols=NCOLS;
nNonZeros=NNONZEROS;
//values=new complex [nNonZeros];
values=new _DOUBLE_COMPLEX_t[nNonZeros];
{for(int i=0; i values.r=real(VALUES);
values.i=imag(VALUES);
cout << " values: i=" << i << ", " << values.r << ", " << values.i << endl;
}}
rowIndex=new int[nRows+1];
{for(int i=0; i rowIndex=ROWINDEX;
}}
columns=new int[nNonZeros];
{for(int i=0; i columns=COLUMNS;
}}
}
SymmSquareMatrix::~SymmSquareMatrix(){
delete [] columns;
delete [] rowIndex;
delete [] values;
}
bool SymmSquareMatrix::Decompose(){
_INTEGER_t error;
MKL_INT opt = MKL_DSS_DEFAULTS;
MKL_INT sym = MKL_DSS_SYMMETRIC_COMPLEX;
MKL_INT type = MKL_DSS_HERMITIAN_INDEFINITE;
{for(int k=0; k cout << "k: " << k << "\t" << rowIndex << endl;
}}
{for(int k=0; k cout << "icol: " << k << "\t" << columns << endl;
}}
//
// Initialize the solver
//
error = dss_create(handle, opt );
if ( error != MKL_DSS_SUCCESS ) {
return false;
}
//
// Define the non-zero structure of the matrix
//
error = dss_define_structure(
handle, sym, reinterpret_cast<_INTEGER_t*>(rowIndex), nRows, nCols,
reinterpret_cast<_INTEGER_t*>(columns), nNonZeros );
if ( error != MKL_DSS_SUCCESS ) {
return false;
}
//
// Reorder the matrix
//
error = dss_reorder(handle, opt, 0);
if ( error != MKL_DSS_SUCCESS ) {
return false;
}
//
// Factor the matrix
//
//error = dss_factor_complex(handle, type, reinterpret_cast<_DOUBLE_COMPLEX_t*>(values) );
error = dss_factor_complex(handle, type, values );
if ( error != MKL_DSS_SUCCESS ) {
return false;
}
return true;
}
bool SymmSquareMatrix::Solve(){
//complex * rhs=new complex [nCols];
//complex * soln=new complex [nCols];
_DOUBLE_COMPLEX_t* rhs=new _DOUBLE_COMPLEX_t[nCols];
_DOUBLE_COMPLEX_t* soln=new _DOUBLE_COMPLEX_t[nCols];
{for(int i=0; i rhs.r=real(RHS);
rhs.i=imag(RHS);
cout << " rhs: i=" << i << ", " << rhs.r << ", " << rhs.i << endl;
}}
_INTEGER_t error;
MKL_INT opt = MKL_DSS_DEFAULTS;
MKL_INT sym = MKL_DSS_SYMMETRIC_COMPLEX;
MKL_INT type = MKL_DSS_HERMITIAN_INDEFINITE;
_INTEGER_t nRhs=1;
//
// Solve
//
//error = dss_solve_complex(handle, opt, reinterpret_cast<_DOUBLE_COMPLEX_t*>(rhs), nRhs, soln );
error = dss_solve_complex(handle, opt, rhs, nRhs, soln );
if ( error != MKL_DSS_SUCCESS ) {
return false;
}
{for(int i=0; i cout << " soln: i = " << i << ", " << soln.r << ", " << soln.i << endl;
}}
delete [] soln;
delete [] rhs;
return true;
}
int main(){
SymmSquareMatrix SymmMat;
bool isOK=false;
isOK=SymmMat.Decompose();
if(!isOK){
return 1;
}
isOK=SymmMat.Solve();
if(!isOK){
return 1;
}
return 0;
}

[Konstantin_A_Intel · Accepted solution]

Hi Vijay,

The problem in this testcase is that you try to solveHermitian matrix, and itmust have real numbers on the diagonal (by the definition of Hermitian matrix), but it does contain non-zero imaginary parts on the diagonal.

For the moment, PARDISO is not able to ignore imaginary part of diagonal elements, so you should explicitly set them to zero.

Sorry for this inconvinience, and I hope this little change ofyour program will not be an issue for you.

Regards,
Konstantin

[vijay1]

Hi Konstantin,

Thanks for the reply. The matrices I solve are typically complex and symmetric (not hermitian, though). So I just changed the MKL DSS options in the earlier code (keeping the rest of the code the same) to the following. Now I get the correct results. Thanks again for your help.

MKL_INT opt = MKL_DSS_DEFAULTS;
MKL_INT sym = MKL_DSS_SYMMETRIC;
MKL_INT type = MKL_DSS_INDEFINITE;

Regards,
Vijay

[Konstantin_A_Intel]

Hi Vijay,

Very good news! Feel free to askhelp in case of any new questions.

Best regards,
Konstantin