Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## PARDISO + ILL CONDITIONED MATRIX

Beginner
1,060 Views

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

4 Replies
Honored Contributor III
1,060 Views
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?
Beginner
1,060 Views
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;

{

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_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_iparm, p->p_msglvl, &ddum, &ddum, &error);

if ( error != 0 )

{

//...

return;

}

*p->p_lLastErrorCode = 0;

}

Honored Contributor III
1,060 Views
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.
Beginner
1,060 Views
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}}