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
- Pardiso error code -4 when using Lagrange multipliers

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

dahlberg

Beginner

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

03-01-2011
10:25 AM

177 Views

Pardiso error code -4 when using Lagrange multipliers

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!

1 Solution

Konstantin_A_Intel

Employee

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

03-01-2011
09:29 PM

177 Views

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

Link Copied

10 Replies

Konstantin_A_Intel

Employee

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

03-01-2011
09:29 PM

178 Views

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

Beginner

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

03-02-2011
02:36 AM

177 Views

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

Employee

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

03-02-2011
09:55 AM

177 Views

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

Employee

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

03-02-2011
10:11 AM

177 Views

Regards,

Konstantin

dahlberg

Beginner

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

03-03-2011
06:41 AM

177 Views

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

Beginner

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

03-08-2011
01:12 AM

177 Views

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

Beginner

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

03-08-2011
01:50 AM

177 Views

** 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

Employee

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

03-08-2011
08:33 PM

177 Views

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

Regards,

Konstantin

dahlberg

Beginner

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

04-29-2011
08:11 AM

177 Views

Hi Konstantin,

Just curious, did you ever get any sense out of that problem?

All the best,

Carl

Just curious, did you ever get any sense out of that problem?

All the best,

Carl

Konstantin_A_Intel

Employee

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

05-02-2011
09:32 PM

177 Views

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

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.