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