Hi!
Your question is not so easy. Despite of PARDISO, DSS doesn't support simultaneous processing of several matrices with one and the same sparse structure. So, it is not possible to examine your situation without working example. Could you provide us with your test case?
If I am not mistaken you want to solve two matrices with one and the same sparse structure at the same time. There are two possible ways for doing this.
First, different matrices with their LU factors can be stored in separate handles. So, you have to reorder, factorize and solve each matrix separately. Of course, in this case memory is used not effectively, because many internal PARDISO arrays are one and the same in both cases.
Second, this problem can be solved by PARDISO which support simultaneous solving of several matrices which have the same sparse structure. This is the best way from the memory utilization point of view. Please see MKL reference manual for details (namely pay attention to maxfct and mnum parameters of PARDISO).
With best regards,
Sergey
Dear Sir,
Thank you for your response and proposal.
My main part of the code for sparse matrix solving is as follows. The dss routines are encapsulated in class CSMat. The function defines sparse matrix object and calls the member function of the object to carry out the sparse matrix calculation. It could be seen in the function that two sparse matrixes, smatA and smatB, are defined at same time. When smatA and smatB are not so large, the function passed. When smatA and smatB are large enough, program crashes at the factorization of the second matrixe.
Is there a possibility to get the LU factorization of the input matrix? So I could store it else where, and reset it when doing solving. Only one sparse matrix would exists in this case.
Thank you very much for your help!
////////////////////////////////////////////
class CSMat
{
public:
CSMat(void);
public:
~CSMat(void);
bool Create8DefineStructure(void);
bool Reorder(void);
bool Factor(void);
bool Solve(int nVec,double* adVec);
bool Solve(CMatrix& m);
bool Delete();
CArray<double,double&> m_adEle;
CArray<int,int&> m_aiRow;
CArray<int,int&> m_aiCol;
private:
_MKL_DSS_HANDLE_t m_handle;
bool m_bDefined;
};
bool CSMat::Create8DefineStructure(void)
{
int nRow,nNonZeros;
int opt = MKL_DSS_DEFAULTS;
int sym = MKL_DSS_SYMMETRIC;
int type = MKL_DSS_POSITIVE_DEFINITE;
int error;
nRow=m_aiRow.GetSize()-1;
nNonZeros=m_aiRow[nRow-1];
error = dss_create(m_handle, opt );
if ( error != MKL_DSS_SUCCESS ) return false;
error = dss_define_structure(
m_handle, sym, m_aiRow.GetData(), nRow, nRow,
m_aiCol.GetData(), nNonZeros );
if ( error != MKL_DSS_SUCCESS ) return false;
m_adEle.SetSize(nNonZeros);
m_bDefined=true;
return true;
}
bool CSMat::Reorder(void)
{
int opt = MKL_DSS_DEFAULTS;
int error;
error = dss_reorder(m_handle, opt, 0);
if ( error != MKL_DSS_SUCCESS ) return false;
return true;
}
bool CSMat::Factor(void)
{
int type = MKL_DSS_POSITIVE_DEFINITE;
int error;
error = dss_factor_real(m_handle, type, m_adEle.GetData());
if ( error == MKL_DSS_OUT_OF_MEMORY ){
AfxMessageBox(_T("OUT_OF_MEMORY "),MB_OK,0);
return false;
}
if ( error != MKL_DSS_SUCCESS ){
AfxMessageBox(_T("ERROR"),MB_OK,0);
return false;
}
return true;
}
bool Function( )
{
double *adGK;
CSMat smatA;
CSMat smatB;
//smatA and smatB coexist …………….
if(!smatA.Create8DefineStructure()){
return false;
}
if(!smatB.Create8DefineStructure()){
return false;
}
//Calculate elements of smatA
adGK=smatA.m_adEle.GetData();
for(loop=0;loop<nEle;loop++){
apEle[loop]->AsmbStf(adGK);
}
//Calculate elements of smatB
adGK=smatB.m_adEle.GetData();
for(loop=0;loop<nEle;loop++){
apEle[loop]->AsmbStf(adGK);
}
if(!smatA.Reorder())
return false;
if(!smatB.Reorder())
return false;
if(!smatA.Factor())
//Ok!!! return false;
if(!smatB.Factor())
//Crashed here !!! return false;
return true;
}
