Does Pardiso do any pivoting?

[Tony_Garratt]

We are giving PARDISO a test driving on some large sparse unsymmetric matrices and using the default values of iparam. However, we are seeing some very large residual norms, which makes me ask the questions (i) does PARDISO do any pivoting and (ii) is there an option we should try to get better numerical results?

[Sergey_K_Intel1]

Quoting - Tony Garratt

We are giving PARDISO a test driving on some large sparse unsymmetric matrices and using the default values of iparam. However, we are seeing some very large residual norms, which makes me ask the questions (i) does PARDISO do any pivoting and (ii) is there an option we should try to get better numerical results?

Dear Tony,

It is not clear from your description whatwasthe value of mtype for calling PARDISO? I assume its value is 11?

In fact, PARDISO uses a pertubed pivot strategy for all types of matrices except mtype=2 (real and symmetric positive), mtype=4 complex and Hermitian positive definite), mtype=1 (real and structurally symmetric) and mtype=4 (complex and structurally symmetric).Iparm(10) instructs PARDISO how to handle small pivots or zero pivots of all other types of matrices. More precisely, iparam(10) is used for computing a threshold. The default values of iparm(10)and more detailed information aboutthe pertubed pivoting strategy can be found in MKL Reference Manual (seethe description of iparm(10)). After factorization, iparm(14) contains the number of pertubed pivots during the elimination process. It is recommended to use iterative refinementfor the solver step (see description of iparm(8)for more details about iterative refinement) if iparm(14) is not zero.Inmany cases the usage of iterative refinement allowsreducing the norm of residual vector.On default, iterative refinement steps is not used.

Large residual norms mean that your matrices areill-conditioned or close to them.I'd recommend toset iparm(11)=1 and iparm(13)=1 for solving such problems.

All the best
Sergey

[Tony_Garratt]

Quoting - Sergey Kuznetsov (Intel)

Dear Tony,

It is not clear from your description whatwasthe value of mtype for calling PARDISO? I assume its value is 11?

In fact, PARDISO uses a pertubed pivot strategy for all types of matrices except mtype=2 (real and symmetric positive), mtype=4 complex and Hermitian positive definite), mtype=1 (real and structurally symmetric) and mtype=4 (complex and structurally symmetric).Iparm(10) instructs PARDISO how to handle small pivots or zero pivots of all other types of matrices. More precisely, iparam(10) is used for computing a threshold. The default values of iparm(10)and more detailed information aboutthe pertubed pivoting strategy can be found in MKL Reference Manual (seethe description of iparm(10)). After factorization, iparm(14) contains the number of pertubed pivots during the elimination process. It is recommended to use iterative refinementfor the solver step (see description of iparm(8)for more details about iterative refinement) if iparm(14) is not zero.Inmany cases the usage of iterative refinement allowsreducing the norm of residual vector.On default, iterative refinement steps is not used.

Large residual norms mean that your matrices areill-conditioned or close to them.I'd recommend toset iparm(11)=1 and iparm(13)=1 for solving such problems.

All the best
Sergey

Thank you very much for your detailed response. mytpe is indeed 11. I will try your suggestions and report back. Im comparing pardiso against some other sparse matrix solvers, so its a good test of pardiso's pivoting and iterative refinement algorithms. Out matrices do have high condition numbers (close to ill-conditioned), but they arise from real engineering applications.

[Tony_Garratt]

Quoting - Sergey Kuznetsov (Intel)

Dear Tony,

It is not clear from your description whatwasthe value of mtype for calling PARDISO? I assume its value is 11?

In fact, PARDISO uses a pertubed pivot strategy for all types of matrices except mtype=2 (real and symmetric positive), mtype=4 complex and Hermitian positive definite), mtype=1 (real and structurally symmetric) and mtype=4 (complex and structurally symmetric).Iparm(10) instructs PARDISO how to handle small pivots or zero pivots of all other types of matrices. More precisely, iparam(10) is used for computing a threshold. The default values of iparm(10)and more detailed information aboutthe pertubed pivoting strategy can be found in MKL Reference Manual (seethe description of iparm(10)). After factorization, iparm(14) contains the number of pertubed pivots during the elimination process. It is recommended to use iterative refinementfor the solver step (see description of iparm(8)for more details about iterative refinement) if iparm(14) is not zero.Inmany cases the usage of iterative refinement allowsreducing the norm of residual vector.On default, iterative refinement steps is not used.

Large residual norms mean that your matrices areill-conditioned or close to them.I'd recommend toset iparm(11)=1 and iparm(13)=1 for solving such problems.

All the best
Sergey

I tried your suggestions and the new iparam(13) (improved accuracy) option in MKL 10.1- it improvesthe residual norm for some of our matrices,but not all. On one in particular casewhere the residual norm is very high, it reducedan order of magnitude, but its still way too large.
For this matrix, I transposed it and solved the system again and the residual norm was much lower and acceptable. Since your backsolve routine does not have an option to solve for the transpose, I get the wrong solution of course. But it does show that there pardiso is having problems with pivoting and/or numerical stability for some of our engineering problems (note that I am comparing pardiso against other third party unsymmetric sparse linear solvers and none of thosehave problems with the residual norm, although they are slower than pardiso generally).

So, I would conclude
(1) There is great value in having a solve transpose iphase option in pardiso
(2) Pardiso may be fine for many matrices, but for a certain class of near to ill-conditioned or ill-conditioned, its pivoting scheme may not beas good as other sparse solvers around

Than

[ArturGuzik]

Quoting - Tony Garratt

(note that I am comparing pardiso against other third party unsymmetric sparse linear solvers and none of thosehave problems with the residual norm, although they are slower than pardiso generally).

(2) Pardiso may be fine for many matrices, but for a certain class of near to ill-conditioned or ill-conditioned, its pivoting scheme may not beas good as other sparse solvers around

Tony,

these other type solvers you're refering to are slower because they implement SVD (as a guard for ill-conditioning) or use exact the same strategy as Pardiso?

A.

[Tony_Garratt]

Quoting - ArturGuzik

Quoting - Tony Garratt

(note that I am comparing pardiso against other third party unsymmetric sparse linear solvers and none of thosehave problems with the residual norm, although they are slower than pardiso generally).

(2) Pardiso may be fine for many matrices, but for a certain class of near to ill-conditioned or ill-conditioned, its pivoting scheme may not beas good as other sparse solvers around

Tony,

these other type solvers you're refering to are slower because they implement SVD (as a guard for ill-conditioning) or use exact the same strategy as Pardiso?

A.

The other two solvers dont use SVD - they are sparse solvers implementing their own version of Guassian Elimination. One is a general purpose nonsymmetric sparse linear solver and the other a frontal method. So they are essentially the same method as Pardiso. However the big difference is how all three solvers are handling fill-in, pivoting and numerical stability. But I amcomparing all solver on matrices arising from a particular class of engineering problems. However, from my testing my conclusion is that Pardiso pivoting strategy is not as good as handling near ill-conditioned problems as others around.

[ArturGuzik]

Tony,

very interesting and really good to know. Many thanks for explanations.

A.

[Sergey_K_Intel1]

Quoting - Tony Garratt

I tried your suggestions and the new iparam(13) (improved accuracy) option in MKL 10.1- it improvesthe residual norm for some of our matrices,but not all. On one in particular casewhere the residual norm is very high, it reducedan order of magnitude, but its still way too large.
For this matrix, I transposed it and solved the system again and the residual norm was much lower and acceptable. Since your backsolve routine does not have an option to solve for the transpose, I get the wrong solution of course. But it does show that there pardiso is having problems with pivoting and/or numerical stability for some of our engineering problems (note that I am comparing pardiso against other third party unsymmetric sparse linear solvers and none of thosehave problems with the residual norm, although they are slower than pardiso generally).

So, I would conclude
(1) There is great value in having a solve transpose iphase option in pardiso
(2) Pardiso may be fine for many matrices, but for a certain class of near to ill-conditioned or ill-conditioned, its pivoting scheme may not beas good as other sparse solvers around

Than

Dear Tony,

According to a detailed analysis of direct sparse solvers done by the CCLRC (UK) (please visit http://epubs.cclrc.ac.uk/bitstream/724/raltr-2005005.pdf), PARDISO pivoting strategy performs very well compared to alternatives.

Ill-conditioned matrices or close to them are quite sensitive to round-off errors, matrix permutations, pivoting strategy and etc. I think the strategy of splitting a matrix into blocks and operating with them are more important for matrices close to ill-conditioned than differences between pivoting strategy implemented in PARDISO and pivoting from other sparse solvers. Other solvers were a little bit lucky in your particular case and they got other internal splittings after reordering. So I'd recommend you to choose another internal PARDISO reordering scheme (the multiple minimum degree method) instead of METIS by setting iparm(2)=1. Internal PARDISO splitting scheme can also be changed by using your own permutation vector and setting iparm(5)=1. The usage of your own permutation vector might help as well. Of course, the ways described above might cause performance degradation and causemore memory consumption. I'd advise you to fix the number of threads in your experiments because the way of splitting matrix into supernodes in PARDISO depends on the number of threads.

Implementation of solve transpose option is in our plans for future MKL releases.

All the best

Sergey

[Tony_Garratt]

Thanks Sergey for the reply. I had completly forgotten about the report that the RAL group did. I have already tried your suggestions (with the exception ofthe permutation - I think its unreasonable of Intel to expect theuser of pardiso to do that, in most cases anyway,since its the job of the sparse solver to do that for them;-), and still I see the high residuals. So either I have a bug in my implementation of how I am calling pardiso or the matrices I am using are "very special".

One question please- when you pass in the matrix do you have to have the column numbers in increasing order (I dont think the manual mentions anything about that)?

If I understand correctly, you are suggesting I should try experimenting with different values of MKL_NUM_THREADS (I have a 2-core machine) to try to get pardiso to break the problem up into more blocks? Are you suggesting also to pass in a permutation (even if it is perm(i) = i, since I dont know what a good one will be apriori) to try to get pardiso to use more blocks?

Would you like me to send you one of our matrices so you can give it a try? It would be embarrassing if I have got a bug in the way I am calling pardiso. If its a bug in my code, I will crawl off into the distance with an embarrassed look on my face and then do my benchmarking again; if there is an issue with pardiso, its a great test case for Intel to work on. Let me know what you think and give me somewhere I can email the matrix to you directly please (I do not want to post it on the forum). One final point - the residual norm I am using is ||b - Ax||/||x||/||A|| (from Golub and Van Loan's "Matrix Compuatations" book).

thanks very much,
Tony

[Sergey_K_Intel1]

Quoting - Tony Garratt

Thanks Sergey for the reply. I had completly forgotten about the report that the RAL group did. I have already tried your suggestions (with the exception ofthe permutation - I think its unreasonable of Intel to expect theuser of pardiso to do that, in most cases anyway,since its the job of the sparse solver to do that for them;-), and still I see the high residuals. So either I have a bug in my implementation of how I am calling pardiso or the matrices I am using are "very special".

One question please- when you pass in the matrix do you have to have the column numbers in increasing order (I dont think the manual mentions anything about that)?

If I understand correctly, you are suggesting I should try experimenting with different values of MKL_NUM_THREADS (I have a 2-core machine) to try to get pardiso to break the problem up into more blocks? Are you suggesting also to pass in a permutation (even if it is perm(i) = i, since I dont know what a good one will be apriori) to try to get pardiso to use more blocks?

Would you like me to send you one of our matrices so you can give it a try? It would be embarrassing if I have got a bug in the way I am calling pardiso. If its a bug in my code, I will crawl off into the distance with an embarrassed look on my face and then do my benchmarking again; if there is an issue with pardiso, its a great test case for Intel to work on. Let me know what you think and give me somewhere I can email the matrix to you directly please (I do not want to post it on the forum). One final point - the residual norm I am using is ||b - Ax||/||x||/||A|| (from Golub and Van Loan's "Matrix Compuatations" book).

thanks very much,
Tony

Dear Tony,

Let me answers to some yourquestions and let's start from the question about column indices. MKL Reference Manual says

ja Integer

Array. ja(*) contains column indices of the sparse matrix A stored in compressed sparse row format. The indices in each row must be sorted in increasing order.

PARDISO in MKL 10.1 and later releases contains new parameter iparm(27) and if iparm(27)=1, the PARDISO also checks whether column indices are sorted in increasing order.

In my previous message I suggested to fix the number of threads while experimenting with PARDISO. I'd propose to set MKL_NUM_THREADS=1.

We have been working onnew PARDISO features for solving ill-conditioned matrices and itwouldbe greattowork withyour matrices. So please visit http://premier.intel.com and submitissue. Please attachyour matrices and a code prototype with PARDISO settings.

Thanks in advance
All the best
Sergey