Pardiso error code -4 when using Lagrange multipliers

[dahlberg]

Hi!

I have built a research FE-code in Fortran 90/95 and I'm using the Pardiso solver for the linear system of equations that results. I have used it for quite some time and it have been working beautifully. The problem I want to solvegives mea real and structurally symmetric matrix, mtype = 1.

The problem is that I recently added the ability to have Lagrange multipliersso I get some zerosalong the diagonal. I get the error code: -4, "zero pivot". I have tried to figure out if I need to supply some additional information to iparm when I have zeros on the diagonal, but to be honest I reachthe end of my linear algebra knowledge once going beyond LU-decomposition...

I have checked that thearrays ia and ja are correct by mimicking the process in matlab and solving the system there - which works just fine, albeit slow and can only handle smaller problems. I set upthe system matrix in matlab using the arrays ia and ja thatI get from the Fortran code, and for small problems I can check and make sure that the resultingmatrixis correct (and I get resonable solutions).

I also tried setting the zeros on the diagonal to a small number, 1.d-30 and 1.d-16 but neither worked.

My current settings to Pardiso are (which I guess should be very close to the default settings):

nrhs = 1

maxfct = 1

mnum = 1

iparm(:) = 0

iparm(1) = 0 ! no solver default

iparm(2) = 2 ! fill-in reordering from METIS

iparm(3) = mkl_get_max_threads() !numbers of processors, value of MKL_NUM_THREADS

iparm(4) = 0 ! no iterative-direct algorithm

iparm(5) = 0 ! no user fill-in reducing permutation

iparm(6) = 0 ! =0 solution on the first n components of x

iparm(7) = 16 ! default logical fortran unit number for output

iparm(8) = 0 ! numbers of iterative refinement steps

iparm(10) = 13 ! perturb the pivot elements with 1E-13

iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS

iparm(14) = 0 ! Output: number of perturbed pivots

iparm(18) = -1 ! Output: number of nonzeros in the factor LU

iparm(19) = -1 ! Output: Mflops for LU factorization

iparm(20) = 0 ! Output: Numbers of CG Iterations

iparm(27) = 1 ! Check matrix formulation (ia, ja and aa)

mtype = 1 ! Matrix type

Any help with this would be greatly appreciated. Many thanks, have a nice day!

[Konstantin_A_Intel · Accepted solution]

Hi,

Did you try to solve this matrix with mtype=11, as unsymmetrical? Anyway, can you attach the testcase to this thread in order we can check if it's a real issue? Or you may just dump (ia, ja, a) arrays in a file in any sensible form.

BTW, Intel MKL PARDISO doesn't support this input parameter:

iparm(7) = 16 ! default logical fortran unit number for output

Another observation is that you use iparm(1)=0 that means you rely on default parameters. If you want to set any parameters manually, please setiparm(1) = 1.

Regards,

Konstantin

[dahlberg]

Hi Konstantin

That helped, thanks alot! I used mtype=11 and now it solves perfectly.

A question relating to performance: in your opinion, should I leave mtype=11 always, or should I set it back to mtype=1 when not using Lagrange multipliers? Is type 11 slower than type 1? Also, is it not a little strange that it did not work for type 1, I mean that is actually the kind of matrix I have and adding Lgr. mult. should not change that since they only add symmetric contributions to the stiffness matrix.

Anyway, thanks again for the quick response.

PS. Yes, the iparm thing was a mistake. I copied that part of the input and I think the previous user might be running the Univ. of Basel version.

Best regards,
Carl

[Konstantin_A_Intel]

Hi Carl,

I'm glad to hear that the workaround worked! Anyway, it would be great to investigate the issue with mtype=1. As far as I'm concerned, you have a specific matrix which PARDISO failed to solve and all was Ok with other matrices, wasn't it? So, could you please provide a reproducer in order I can investigate the root cause of the failure?

Re perfromance aspect.. indeed, algorithms for 1 and 11 types a little bit different, so the perfromance may differ to some extent. But I beleive the difference won't be outstanding and recommend using mtype=11 for a while.

Regards,

Konstantin

[Konstantin_A_Intel]

Another thought.. you said that the matrix has a lot of zeros on the diagonal - this can be a problem for mtype=1. First of all, this type has strong requirement for input matrix - namely, no implicit zeros are allowed on the main diagonal (in other words, you must set all the diagonal elements as explicit zero values in the matrix). Unsymmetrical type doesn't have such a limitation. And one more thing - unsymmetrical type has more comprehensive pivoting strategy that may impact the stability.

Regards,

Konstantin

[dahlberg]

Hi Konstantin

I will try to get you a sample matrix as soon as possible. Code is working now so the focus is on getting results right now... I will put it high up on the "to do list".

You did understand correct. I had a matrix, lets call it K, that would solve with mtype = 1. Then I augmented that with some Lagrange multipliers given by the matrix C (almost all zeros and a few 1 and -1)so that the new system is

K_aug = [ K C' ]
[ C 0 ]

where C' denotes the transpose of C, and 0 should be interpreted as the zero-matrix of apropriate dimensions. I made sure that the zeros on the diagonal was explicitly stored.

I did quite a few test-runs last night and in most cases the solver worked perfectly once I switched to mtype = 11. However I still occationally get the "zero pivot" error, but only in extreme cases and I can skip around that by reducing my step size in the equilibrium iterations and playing around a little with iparm(4) and iparm(21).

Anyway, once again thanks for the help. I'll get back to you with some test cases.

Regards,
Carl

[dahlberg]

Hi Konstantin

I made two example structures, one with and one without Lagrange multipliers. One that solves with mtype = 1 and one that won't solve unless I change to mtype = 11. Except for the added Lagrange multpliers they are identical.

I tried to add them via the "Add Files" button, but I must say that the interface was a little confusing. I only got to create a folder, and I couldn't see any way of actually adding the files. How do I do this? Would you prefer to get them via mail instead?

Regards,
Carl

[dahlberg]

Ok, think I figured out how to do it now. Anyway, the files are pure text-files. The structure is; two stars **marks the start of a new array, followed by the name (using standard naming), and then a number to denot the length of the array. For example:
** ia, 206
29
...
** ja, 9313
1
2
...
and so forth. I supplied ia, ja, aa, and rhs.

Sorry for the slightly retarded format, but it felt like the simplest method. The arrays correspond to a fairly small problem.

The array aa might turn out to be symmetric (i.e. mtype = 2), but that is only because it comes from the first iteration of the algorithm. Later on in the solution process the stiffness matrix is guaranteed to be only strucurally symmetric (mtype = 1).

The structure in Arrays_Standard.out solves with mtype = 1, but in Arrays_LgrMult.out the structureis augmented with 16 Lagrange multiplier equations, and now it only solves with mtype = 11. But the type of matrixhave not have changed. Zeros on the diagonal are explicitly stored.

Regards,
Carl

[Konstantin_A_Intel]

Carl, thank you! I''ll grab the files.

Regards,

Konstantin

[dahlberg]

Hi Konstantin,
Just curious, did you ever get any sense out of that problem?
All the best,
Carl

[Konstantin_A_Intel]

Hello Carl,

I have reproduced your test, but unfortunately did not have a time to deeply investigate it. But as I mentioned in some of the previous posts, the behavior is quite expected for mtype=1. Look, this type doesn't support many techniques available for mtype=11 which is needed to enhance reliability of the numerical factorization process(such as matching and scaling). So, I would not recommend you to use mtype=1 at all, especially for highly indefinite systems. Please use mtype=11.

Regards,

Konstantin