PARDISO: different outputs from the same source

7nic9 · ‎11-10-2008

Dear all,

I'm using PARDISO to solve a linear system originated from a Stokes problem. The matrix is real and structurally symmetric.

My first problem is that if I treat it as structurally symmetric (mtype=1) the solver fails and returns:

*** error PARDISO ( numerical_factorization) error_num= -1
*** error pardiso: zero pivot

.........

PARDISO statistics

......

error in PARDISO = -4

However, if I treat it as unsymmetric (mtype=11) the solver does not return any error. But perhaps this is not even a problem and depends on the fact that two different methods are used for these two type of matrices.

The second (and more serious) problem is that, using mtype=11, running the same executable twice results in two different solutions. To give an idea, for an element of the solution vector of the order of unity, the differences between the two runs can be up to 1e-5 in double precision. I use standard options (iparm(1)=0). I also tried to use customized options, but (strangely) no matter if I use iparm(2)=0 or iparm(2)=2, after the reordering phase (phase=11) Pardiso returns what follow:

*** error PARDISO (prereordering_mtype11,13) error_num= 0
*** error pardiso: preordering failed after 0 neqns out of 17312
struct. sing. or input/parameter problem (matrixtype 11,13)

.........

PARDISO statistics

......

error in PARDISO = -6

Does anybody have any idea of what can be wrong?

Thank you very much in advance,

Nicola

Sergey_K_Intel1 · ‎11-13-2008

Dear Nicola,

You are right, two different methods are used for mtype=1 and mtype=11.

The settings iparm(2)=0 and iparm(2)=2 use the same type of reordering namely METIS reordering. There exists another built-in reordering scheme so called MMD reordering available through iparm(2)=1. So you could try to useiparm(2)=1 instead of METIS.

There exists the following limitation in PARDISO usage (see the description of array ja, pp 2352 in MKL 10.1 Gold Reference Guide):

Ja - 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.

I suspect PARDISO fails in your casebecause the representation of a matrix in the CSR format is wrong.

MKL 10.1 Gold provides the new feature for checking whether this condition is satisfied. The checker is turn off on default. To turn it on, iparm(27) must be set to 1. Here is what MKL reference manual says

If iparm(27)=1, then PARDISO check integer arrays ia and ja. In particular, PARDISO checks whether column indices are sorted in increasing order within each row.

All the best

Sergey

7nic9 · ‎11-18-2008

Dear Sergey,

thank you very much for your answer. I think the arrays are stored correctly. I just tried the new MKL 10.1 with iparm(27)=1 and I do not receive any error. Actually, I finally managed to get always the same output using the unsymmetric solver (mtype=11) with the following options:

iparm(1)=1
iparm(2)=2
iparm(3)=mkl_get_max_threads()
iparm(4)=0
iparm(5)=0
iparm(6)=0
iparm(8)=-5
iparm(10)=7
iparm(11)=1
iparm(13)=1
iparm(18)=-1
iparm(19)=0
iparm(20)=0
iparm(27)=1

I realized that the key problem was the value of iparm(10). By trial and error, I found that perturbing the pivot with values smaller than 10^-7 results in a not very precise solution that can change from run to run. With these options my code works fine and I'm happy with it. However, although the matrix is structurally symmetric, using mtype=1, the solver still fails to factorize the matrix. After phase=11 and, apparently, no matter what options I use (even using iparm(2)=1) this is what I get for a typical run

================ PARDISO: solving a real struct. sym. system ================

Summary PARDISO: ( reorder to reorder )
================

Times:
======

Time fulladj: 0.001180 s
Time reorder: 0.085897 s
Time symbfct: 0.009929 s
Time malloc : -0.002462 s
Time total : 0.111708 s total - sum: 0.017165 s

Statistics:
===========
< Parallel Direct Factorization with #processors: > 1
< Numerical Factorization with Level-3 BLAS performance >

< Linear system Ax = b>
#equations: 10858
#non-zeros in A: 96286
non-zeros in A (%): 0.081670
#right-hand sides: 1

< Factors L and U >
#columns for each panel: 128
#independent subgraphs: 0
< Preprocessing with state of the art partitioning metis>
#supernodes: 7605
size of largest supernode: 176
number of nonzeros in L 298753
number of nonzeros in U 231369
number of nonzeros in L+U 530122
*** error PARDISO ( numerical_factorization) error_num= -1
*** error pardiso: zero pivot

With best regards,

Nicola

Gennady_F_Intel · ‎11-26-2008

I would recommend you submit the issue against MKL to Premier support( https://premier.intel.com/ ). I will work with this issue there.--Gennady

Sergey_K_Intel1 · ‎11-26-2008

Dear Nicola,

Since PARDISO found zero pivots, your system is ill-conditioned or singular. In this case, I'd suggest to experiment with different regularization strategies like Tikhonov regularization (see http://en.wikipedia.org/wiki/Ridge_regression ) or any other before calling PARDISO. The purpose of regularization is to introduce additional information in order to solve an ill-posed problem.

All the best

Sergey

7nic9 · ‎11-28-2008

Quoting - Gennady Fedorov (Intel)

Thank you. I submitted the issue to the Premier Support.

Nicola

I would recommend you submit the issue against MKL to Premier support( https://premier.intel.com/ ). I will work with this issue there.--Gennady

7nic9 · ‎11-28-2008

Quoting - Sergey Kuznetsov (Intel)

Dear Sergey,

thanks for the hint. I'll try to have a look at the Tikhonov regularization and see if I manage to improve my matrix conditionality. By the way, dose Pardiso (or some other mkl tool) allow to estimate the condition number?

Best regards,

Nicola

Dear Nicola,

Since PARDISO found zero pivots, your system is ill-conditioned or singular. In this case, I'd suggest to experiment with different regularization strategies like Tikhonov regularization (see http://en.wikipedia.org/wiki/Ridge_regression ) or any other before calling PARDISO. The purpose of regularization is to introduce additional information in order to solve an ill-posed problem.

All the best

Sergey