Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Problem with Pardiso and indefinite matrix

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Moritz_K_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-13-2017
04:34 AM

383 Views

Problem with Pardiso and indefinite matrix

Hello,

I am using Pardiso to solve linear equations with an symmetric indefinite matrix (saddle-point structure). Usually this works quite good, however sometimes the solver teminates with the error "*** error PARDISO: iterative refinement contraction rate is greater than 0.9, interrupt".

Attached you find the csr representation of my matrix (symmetric indefinite so only the upper right triangle is represented). Matlab tells me that the condition of the matrix is 4.596e+05 and is able to solve the equation. Below you find the settings I used for the Pardiso call. In particular I am using METIS.

After some searching I stumbled upon a post with a similar problem where it was proposed to try a different reordering strategy for badly conditioned matrices. Switching to the minimum degree algorithm helped for this matrix. However I find it strange that METIS has a problem with this matrix since the condition is not too bad. Are the rest of my parameter choices sensible?

Help is greatly appreciated!

Best Moritz

MKL_INT mtype = -2;

MKL_INT nrhs = 1;

void *pt[64];

/* Pardiso control parameters. */

MKL_INT iparm[64];

MKL_INT maxfct, mnum, phase, error, msglvl;

char *uplo;

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

/* .. Setup Pardiso control parameters. */

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

for (i = 0; i < 64; i++)

{

iparm* = 0;
}
iparm[0] = 1;
iparm[1] = 2;
iparm[3] = 0;
iparm[4] = 0;
iparm[5] = 0;
iparm[6] = 0;
iparm[7] = 10;
iparm[8] = 0;
iparm[9] = 8;
iparm[10] = 1;
iparm[11] = 0;
iparm[12] = 1;
iparm[13] = 0;
iparm[14] = 0;
iparm[15] = 0;
iparm[16] = 0;
iparm[17] = 0;
iparm[18] = 0;
iparm[19] = 0;
iparm[20] = 1;*

iparm[24] = 1;

iparm[26] = 1;

iparm[34] = 0;

maxfct = 1;

mnum = 1;

msglvl = 1;

error = 0;

uplo = "non-transposed";

Link Copied

10 Replies

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-13-2017
09:21 AM

383 Views

Hi,

That's well knows problem that caused by complexity of full pivoting strategy. In few words if you want implement LU decompositon for following matrix:

| 0 1 |

|1 1 |

you will got division by zero despite good condition number. For huge real matrix the situation is slightly different and "zero" pivots replace on eps during factorization but this trick can significantly broke resulted LU matrix and caused situation that you have

Thanks,

Alex

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-13-2017
09:27 AM

383 Views

Moritz_K_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-13-2017
09:54 AM

383 Views

@Alex: perhaps I do not understand something fully, but we do we not use a LDL-factorization for a symmetric indefinite matrix? Can you comment on the rest of my parameter choices for this kind of problem? Thanks!

@mecej4: The matrix is symmetric (indefinite), so in accordance with the requirements of pardiso only the upper triangle is represented in the csr-format. This upper triangle contains zero columns/ rows, however the full matrix has full rank.

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-13-2017
10:46 AM

383 Views

Hi,

The same situation for LDL^t and LU decomposition, i've described LU only as example.

Parameters looks fine, but what is the size of your problem? I am bit confused why did you set sequential solving step?

Thanks,

Alex

Moritz_K_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-14-2017
05:27 AM

383 Views

Hey Alex,

ok thanks! The sequential solving step is a relict of trying to eliminate possible error reasons.

While I am now satisfied with the solution for the original matrix (relative residuum 1e-14) the situation for a model with finer discretization leading to the attached matrix is much worse. Calculating the relative residuum of Ax-b gives 2.161726e-008 which is quite large for a direct solver...

Best Moritz

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-14-2017
05:41 AM

383 Views

Hi Moritz,

Thanks for the matrix, I will test it but if you have time you can also implement following trick on your side. Let's turn on Schur complement calculation

1.set iparm[35] to 2

2.fill perm to 0 if this position doesn't corespond to zero diagonal element and to 1 otherwise

3.allocate array x to max(zero_block*zero_block, neqns)

4. call pardiso with separate calls, ido=11,22,33,-1

From my perspective in such solution vector will be quite good

Thanks,

Alex

philippedevloo

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-08-2020
07:11 AM

346 Views

I am having the same problem for inverting matrices originated from saddle point problems. I understand your suggestion was to "move" all Lagrange multipliers to the end of the linear system and use "static condensation" to "safely" invert the matrix.

I would like to suggest a second strategy, but need your expertise to implement it. By slightly changing the permutation vector suggested by pardiso, I can "guarantee" the system can be decomposed without encountering a zero pivot. Therefore, I would need to "get" the permutation vector, modify it and "feed it back" to pardiso.

Can you help?

Kirill_V_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-08-2020
09:41 AM

334 Views

I am afraid Alex will no longer help you here.

I see your idea about modifying the permutation vector. So, you can call pardiso phase 11 and get the perm back by setting iparm[4]=1. Then you can modify it as you see fit and call pardiso (better with a new handle, I would suggest destroying the first handle after you get the permutation back), phase 13 e.g.

Where do you see difficulties in doing it this way? Or maybe I misunderstood your idea, feel free to correct me.

Best,

Kirill

philippedevloo

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-10-2020
07:40 AM

320 Views

@Kirill_V_Intel My doubt is in reference to this observation from the users guide

"Setting iparm[4]

Does that mean I will lose "most of the efficiency" of the solver if I supply the permutation vector?

Thank you for your reply!!

Kirill_V_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-17-2020
03:35 PM

286 Views

Could you provide your matrix type and iparm settings? I'd like to double check the cited statement from the documentation. It might have been a bit too cautios.

The idea behind that statement is that much of the parallelism comes from the fill-in reducing reordering and if you supply your permutation, we cannot guarantee that it will work well since you don't allow PARDISO do its own reordering .

However, unless we have a blocker somewhere (this is what I want to check for your iparm settings), if you take a permutation from PARDISO and modify it the way which does not affect much the fill-in and dependency tree, performance should remain the same.

Best,

Kirill

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

For more complete information about compiler optimizations, see our Optimization Notice.