- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Thanks,
Vijay
1 Solution
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
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
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
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
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
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
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;
}
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
values.i=imag(VALUES);
cout << " values: i=" << i << ", " << values.r << ", " << values.i << endl;
}}
rowIndex=new int[nRows+1];
{for(int i=0; i
}}
columns=new int[nNonZeros];
{for(int i=0; i
}}
}
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
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
}}
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
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
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
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
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
Hi Vijay,
Very good news! Feel free to askhelp in case of any new questions.
Best regards,
Konstantin
Very good news! Feel free to askhelp in case of any new questions.
Best regards,
Konstantin
Reply
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