- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Thanks,

Vijay

1 Solution

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Link Copied

6 Replies

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

//double RHS[NCOLS] = { 1, 2, 3, 4, 5 };

complex

class SymmSquareMatrix{

public:

int nRows;

int nCols;

int nNonZeros;

//complex

_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

values=new _DOUBLE_COMPLEX_t[nNonZeros];

{for(int i=0; i

*.r=real(VALUES*

*);*

valuesvalues

*.i=imag(VALUES**);*

cout << " values: i=" << i << ", " << valuescout << " values: i=" << i << ", " << values

*.r << ", " << values**.i << endl;*

}}

rowIndex=new int[nRows+1];

{for(int i=0; i rowIndex}}

rowIndex=new int[nRows+1];

{for(int i=0; i

*=ROWINDEX**;*

}}

columns=new int[nNonZeros];

{for(int i=0; i columns}}

columns=new int[nNonZeros];

{for(int i=0; i

*=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 }}

}

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

}}

{for(int k=0; k

}}

//

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

//complex

_DOUBLE_COMPLEX_t* rhs=new _DOUBLE_COMPLEX_t[nCols];

_DOUBLE_COMPLEX_t* soln=new _DOUBLE_COMPLEX_t[nCols];

{for(int i=0; i

*.r=real(RHS*

*);*

rhsrhs

*.i=imag(RHS**);*

cout << " rhs: i=" << i << ", " << rhscout << " 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}}

_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

*.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;

}

}}

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;

}

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Best regards,

Konstantin

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page