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

Problem with Pardiso and indefinite matrix

Moritz_K_
Beginner
1,241 Views

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

 

 

 

0 Kudos
10 Replies
Alexander_K_Intel2
1,241 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

0 Kudos
mecej4
Honored Contributor III
1,241 Views

Unless I am making a mistake, I see a major problem with the matrix data. The matrix is upper triangular, and rows 49 to 72 are zero-filled. That makes the matrix singular.

0 Kudos
Moritz_K_
Beginner
1,241 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.

0 Kudos
Alexander_K_Intel2
1,241 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 

0 Kudos
Moritz_K_
Beginner
1,241 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

0 Kudos
Alexander_K_Intel2
1,241 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

0 Kudos
philippedevloo
Beginner
1,204 Views

@Alexander_K_Intel2 

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?

0 Kudos
Kirill_V_Intel
Employee
1,192 Views

Hi @philippedevloo 

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

0 Kudos
philippedevloo
Beginner
1,178 Views

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

"Setting iparm[4] = 1 prevents use of a parallel algorithm for the solve step."

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

Thank you for your reply!!

 

0 Kudos
Kirill_V_Intel
Employee
1,144 Views

@philippedevloo,

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

0 Kudos
Reply