PARDISO + ILL CONDITIONED MATRIX

Alemdar__Bulent · ‎03-05-2012

Hi All

I am using PARDISO to solve a system of KX=F (K: real symmetric definite matrix).

I have these two cases: K1 and K2 (both matrices are given below). PARDISO returns a success for K1 and it gives an error for K2 (instability in the system).

In fact, both K1 and K2 are not well-conditioned matrices (i.e., they are ill-conditioned) but interestingly the PARDISO returns success for K1. The condition numbers I calculated are something like 10^16 for K1 and 10^17 for K2.

Is there a way to avoid this kind of problems? Any way to set some parameters for PARDISO that catches this? We have other (in-house) solvers that catches this, and I am suprised that PARDISO fails to detect it !..

Regards

Bulent

P.S. My system configuration: Windows7, 64bit, VS2010 C++

K1:

%%MatrixMarket matrix coordinate real symmetric

8 8 16

1 1 14098.547767242315000000

1 3 3485.242905845857900000

1 5 6812.909667832833300000

2 2 14098.547767242315000000

2 6 3485.242905845857900000

2 8 6812.909667832833300000

3 3 82377.582791945460000000

3 5 3485.242905845857900000

3 6 -81215.835156663510000000

4 4 81215.835156663510000000

4 7 -0.000000000000113687

5 5 14098.547767242315000000

6 6 82377.582791945460000000

6 8 3485.242905845857900000

7 7 81215.835156663510000000

8 8 14098.547767242315000000

K2:

%%MatrixMarket matrix coordinate real symmetric

8 8 17

1 1 14571.276198818963000000

1 3 3642.819049704741200000

1 5 7285.638099409481600000

2 2 14571.276198818963000000

2 6 3642.819049704741200000

2 8 7285.638099409481600000

3 3 82430.108173231754000000

3 5 3642.819049704741200000

3 6 -81215.835156663510000000

4 4 81215.835156663510000000

4 5 -0.000000000000454747

4 7 0.000000000000113687

5 5 14571.276198818963000000

6 6 82430.108173231754000000

6 8 3642.819049704741200000

7 7 81215.835156663510000000

8 8 14571.276198818963000000

mecej4 · ‎03-05-2012

Trying your matrices out using Matlab's sparse matrix, (and taking what you gave as the upper triangle of a symmetric matrix), I found the estimated condition numbers of both matrices to be about 10.

Something is inconsistent in the story, or I have made a mistake. Would you please post the code that uses Pardiso on these matrices?

Alemdar__Bulent · ‎03-05-2012

Thanks for the quick reply. That is correct that only the upper diagonal terms are given, and the matrices are symmetric. I used Mathematica to calculate condition numbers but I might have done something wrong for calculation of condition numbers. Nevertheless, my routine returns a success with K1 (we have been using INTEL-PARDISO last several years and this is the first time I run into this).

Anyway, below is my code that uses PARDISO.

Regards
Bulent

m_iparm[0] = 1;

m_iparm[1] = 3;

m_iparm[2] = 0;

m_iparm[3] = 0;

m_iparm[4] = 0;

m_iparm[5] = 1;

m_iparm[7] = 2;

m_iparm[8] = 0;

m_iparm[9] = 13;

m_iparm[10] = 1;

m_iparm[11] = 0;

m_iparm[12] = 0;

m_iparm[13] = 0;

m_iparm[14] = 0;

m_iparm[15] = 0;

m_iparm[16] = 0;

m_iparm[17] = -1;

m_iparm[18] = -1;

m_iparm[19] = 0;

m_iparm[20] = 0;

m_iparm[24] = 1;

m_iparm[26] = 0;

m_iparm[27] = 0;

m_iparm[59] = 0;

void ThreadedSolver(void* param)

{

double ddum; /* Double dummy*/

int idum; /* Integer dummy.*/

int* Perm = NULL;

PARDISO ( (_MKL_DSS_HANDLE_t*) p->p_pt, p->p_maxfct, p->p_mnum, p->p_mtype, &phase,

p->p_NumberOfEqns,

p->p_GStiffnessVector, p->p_RowIndexVector, p->p_ColumnsVector,

&idum,

p->p_lNumLoadVectors,

p->p_iparm,

p->p_msglvl,

&ddum, &ddum, &error);

if ( Perm != NULL )

delete [] Perm;

if (error != 0)

{

//..

return;

}

// --------------------------------------------------------------------*/

// .. Numerical factorization.*/

// --------------------------------------------------------------------*/

phase = 22;

PARDISO ( (_MKL_DSS_HANDLE_t*) p->p_pt, p->p_maxfct, p->p_mnum, p->p_mtype, &phase,

p->p_NumberOfEqns,

p->p_GStiffnessVector, p->p_RowIndexVector, p->p_ColumnsVector,

&idum,

p->p_lNumLoadVectors,

p->p_iparm, p->p_msglvl, &ddum, &ddum, &error);

if ( error != 0 )

{

//...

return;

}

*p->p_lLastErrorCode = 0;

}

mecej4 · ‎03-05-2012

I modified the example pardiso_sym_c.c included with MKL 10.3.9 to run your examples, with all right hand side elements b = 1. Both examples ran to completion with no error messages.

I repeated the runs with the Pardiso 4.1.2 (not from Intel) and again both examples ran without errors.

Alemdar__Bulent · ‎03-12-2012

Hi Mecej4
Thanks for having time to look into this. I am a bit suprised with your results. When I run these matrices in my code with MKL-PARDISO, I am getting success for K1 and failure for K2. In fact, I was expecting no success for K1 too.

These matrices are from a finite element model that is intentionally unstable and hence, it was expected that PARDISO should have reported numerical instability for both K1 and K2.

I just run the matrix "K1" in Mathematica (i.e., taking inverse of K1). It returns some results but it also indicates that "Inverse::luc: Result for Inverse of badly conditioned matrix", which was expected (see the log below)

Please note that I am using matrix type "2" for PARDISO solution (i.e., real symmetric definite matrix).

Regards
BUlent

Mathematica Results:

K1 = {{14098.547767242315, 0, 3485.242905845858, 0, 6812.909667832833, 0, 0, 0}, {0, 14098.547767242315, 0, 0, 0, 3485.242905845858, 0, 6812.909667832833}, {3485.242905845858,

0, 82377.58279194546, 0, 3485.242905845858, -81215.83515666351, 0, 0}, {0, 0, 0, 81215.83515666351, 0, 0, -1.13687*10^-13, 0}, {6812.909667832833, 0, 3485.242905845858, 0,

14098.547767242315, 0, 0, 0}, {0, 3485.242905845858, -81215.83515666351, 0, 0, 82377.58279194546,

0, 3485.242905845858}, {0, 0, 0, -1.13687*10^-13, 0, 0, 81215.83515666351, 0}, {0, 6812.909667832833, 0, 0, 0, 3485.242905845858, 0, 14098.547767242315}};

Inverse[K1]

Inverse::luc: Result for Inverse of badly conditioned matrix {{14098.5,0.,3485.24,0.,6812.91,0.,0.,0.},{0.,14098.5,0.,0.,0.,3485.24,0.,6812.91},{3485.24,0.,82377.6,0.,3485.24,-81215.8,0.,0.},{0.,0.,0.,81215.8,0.,0.,-1.13687*10^-13,0.},{6812.91,0.,3485.24,0.,14098.5,0.,0.,0.},{0.,3485.24,-81215.8,0.,0.,82377.6,0.,3485.24},{0.,0.,0.,-1.13687*10^-13,0.,0.,81215.8,0.},{0.,6812.91,0.,0.,0.,3485.24,0.,14098.5}} may contain significant numerical errors. >>

{{5.55309*10^9, 5.55309*10^9, -3.33185*10^10, 0.,

5.55309*10^9, -3.33185*10^10, 0., 5.55309*10^9}, {5.55309*10^9,

5.55309*10^9, -3.33185*10^10, 0., 5.55309*10^9, -3.33185*10^10, 0.,

5.55309*10^9}, {-3.33185*10^10, -3.33185*10^10, 1.99911*10^11,

0., -3.33185*10^10, 1.99911*10^11, 0., -3.33185*10^10}, {0., 0., 0.,

0.0000123129, 0., 0., 1.72357*10^-23, 0.}, {5.55309*10^9,

5.55309*10^9, -3.33185*10^10, 0., 5.55309*10^9, -3.33185*10^10, 0.,

5.55309*10^9}, {-3.33185*10^10, -3.33185*10^10, 1.99911*10^11,

0., -3.33185*10^10, 1.99911*10^11, 0., -3.33185*10^10}, {0., 0., 0.,

1.72357*10^-23, 0., 0., 0.0000123129, 0.}, {5.55309*10^9,

5.55309*10^9, -3.33185*10^10, 0., 5.55309*10^9, -3.33185*10^10, 0.,

5.55309*10^9}}