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

Problem with PARDISO

Danesh_Daroui
Beginner
2,118 Views
Hi,

I am trying to use PARDISO (in Intel C++ compiler 12.0) as sparse direct solver in our application. For testing purposes, I am first assembling the coefficient matrix in the normal dense form and solve it using normal direct solver (getrf, getrs) then I convert it into CSR format and again solve the equation using PARDISO to get the answer and compare it with two solvers. PARDISO results match very well with direct solvers, but in some problems they are totally different and I know that PARDISO returns wrong results. I save the coefficient matrix "M" in the equation which is non-symmetric and in complex format and is a "Msize x Msize" matrix:

Mx=b

and then convert it to CSR format as below:

int nnz=0;

for (int i=0; i if (M!=(dcmx)0.0)
++nnz;

dcmx *a;
int *ia;
int *ja;

a=new dcmx[nnz];
ia=new int[Msize+1];
ja=new int[nnz];

int p1=0;
int p2=0;

for (int i=0; i{
bool f=true;

for (int j=0; j {
if (M[j*Msize+i]!=(dcmx)0.0)
{
// first element in this row is detected
if (f)
{
ia[p2]=p1+1;
++p2;
f=false;
}

a[p1]=M[j*Msize+i];
ja[p1]=j+1;
++p1;
}
}
}

ia[p2]=p1+1;



when formatting is done, I call PARDISO with the following parameters:

iparm[0]=1; // no solver default
iparm[1]=2; // fill-in reordering from parallel METIS
iparm[2]=1; // not used
iparm[3]=0; // use direct solver
iparm[4]=0; // no user fill-in reducing permutation
iparm[5]=0; // Write solution into x
iparm[6]=0; // not used
iparm[7]=2; // max numbers of iterative refinement steps
iparm[8]=0; // not used
iparm[9]=13; // Perturb the pivot elements with 1E-13
iparm[10]=1; // Use nonsymmetric permutation and scaling MPS
iparm[11]=0; // not used
iparm[12]=1; // maximum weighted matching algorithm is switched-on (default for non-symmetric)
iparm[13]=0; // Output: Number of perturbed pivots
iparm[14]=0; // not used
iparm[15]=0; // not used
iparm[16]=0; // not used
iparm[17]=-1;// output: number of nonzeros in the factor LU
iparm[18]=-1;// output: Mflops for LU factorization
iparm[19]=0; // output: numbers of CG Iterations

maxfct=1; // maximum number of numerical factorizations
mnum=1; // which factorization to use
msglvl=1; err=0;


As I stated some times PARDISO results wrong solution and when this happens and I re-run the tests I get different answers from PARDISO, even different from before and that's very odd! Can you please help me to verify where the formatting to CSR or PARDISO parameters that I use are incorrect or there can be something else?

Regards,

D.
0 Kudos
1 Solution
Konstantin_A_Intel
2,118 Views
Hi Danesh,

It seems it was an error in my computations in the previous post - I passed mistakenly a bit changed matrix to SVD routine. I'm sorry! In fact, there're no zero singular values in the matrix and hence it's NOT singular! Condition number is quite big (~10^14). However, it should be enough to produce more or less sensible solution.
We'll look closely into the problem and keep you informed when any results of investigation will be ready.
Thank you very much for this report!
Regards,
Konstantin

View solution in original post

0 Kudos
19 Replies
Konstantin_A_Intel
2,118 Views
Hi Danesh,
The most possible thing is that some of your matrices are ill-conditioned.. that would make the process of solving unstable. What is the size of your matrix? If it's not very large, you can provide us with the matrix and we would estimate its condition number.
Regards,
Konstantin
0 Kudos
Konstantin_A_Intel
2,118 Views
Another note...
As you mentioned you're to use LAPACK, you can call SVD for your matrix and compute conditional number:
\kappa(A) = \frac{\sigma_\max(A)}{\sigma_\min(A)}wheremax(A)andmin(A)are maximal and minimalsingular valuesofArespectively,
Regards,
Konstantin
0 Kudos
Danesh_Daroui
Beginner
2,118 Views
Hi and thanks for your answer. Well, I know that my matrices are highly ill-conditioned, but the questions is why should PARDISO bother? I mean, PARDISO is a direct solver and ill-conditioned matrices should not affect its performance. When iterative solvers are used, then yes, a preconditioner should be used, but does it really matter when a sparse direct solver like PARDISO is used? I even didn't use PARDISO's built-in CG solver. The most tedious task in iterative solvers is to find an appropriate preconditioner! If PARDISO may fail with an ill-conditioned solver, so why do we use PARDISO and not GMRES or CG? :)

On the other hand, as I stated, I get different results when I run same problem! That is very very strange! It feels like the solver is not thread safe or something like that. I will try to post my input matrix and a sample code here so you can run your self and see where the problem is.

Regards,

D.
0 Kudos
Konstantin_A_Intel
2,118 Views
Hi Danesh,
When I asked you either the matrix is ill-conditioned or not I meant that exactly stability (not performance) may hurt in this case. Of course, PARDISO direct solver will produce some results without performance degradation for the matrix with a condition number, say, more than 10^16 - but there's no any guarantee that a solution is precisely the "solution".
You can find rather interesting and simple paper about condition number here:
As far as I'm concerned about iterative solvers, it seems that the task of finding a good preconditioner is not a simple task.
Regards,
Konstantin
0 Kudos
Danesh_Daroui
Beginner
2,118 Views
Hi,

As I said, I expect that my coefficient matrix to be ill-conditioned, but the question is that why PARDISO should fail even if my matrix is highly ill-conditioned? PARDISO is a direct solver! It factorizes the matrix and solves it using backward/forward substitution and that's all! I know that pivoting can be (very very rare) the source of some problems in massively parallel codes which is not an issue here and in that case diagonal scaling may help.

Anyway, I can *not* understand why PARDISO should fail if my system is ill-conditioned! Direct solvers are supposed to be independent from condition number of the system because the solution will not converge like iterative solvers that you can guess number of iterations from spectrum on the eigenvalues. It solves the problem in one shot like all other direct solvers e.g. getrf, getrs in LAPACK! If not, why it is called direct solver? If my system is singular and therefore not solvable, so why LAPACK can solve it?

The other issue is, why the solution differs each time I solve my system with exactly same input data? That is very odd!

I will post my input matrix with the source code so you can experiment it yourself.

D.
0 Kudos
Sergey_K_Intel1
Employee
2,118 Views

Hi Danesh

It is also well known (I highly recommend to read the article by David Goldberg "What Every Computer Scientist Should Know About Floating-Point Arithmetic") that if we change an algorithm in a program with floating point arithmetic, the results of the program will also change. It should be noted that PARDISO uses dynamic parallelization technique. That means that the order of operations is determined dynamically at run time and as a consequence the order of operations might varies from run to run. That leads to the deviation in the results since round-off errors might be different.

The condition number is, in fact, the measure of sensitivity of the solution to round-off errors. There exists estimates for the admissible errors in solution of linear systems. These estimates are well known and can be found in books on linear system solvers (see for example the book G.W.Stewart Matrix Computations") and these estimates are based on a careful analysis of impact of round-off errors on the solution. So lets' take a simplest estimate which is the following

|| (exact solution) - (computed solution) || <= eps*( the condition number of the system) * || right hand side||

where eps is the relative precision for a given type of input data. For example in the case of single complex arithmetic operation eps is 1.1920929E-07, for double precision eps is about 1^(-16).So the larger condition number is, the larger deviations might be.

All the best
Sergey

0 Kudos
Danesh_Daroui
Beginner
2,118 Views

I see. So is there any solution for my case? The interesting thing is that I get different results when the solution, given by PARDISO is already wrong! In all other cases, PARDISO results same answer as many times as I run it and I have verified those answers and I am sure that they are correct.

Besides, how can high condition number affect the results in PARDISO while LAPACK solvers return correct results? Anyway, can there be any solution for my case?

Regards,

D.
0 Kudos
Konstantin_A_Intel
2,118 Views
Danesh,
It's a possible situation when LAPACK is able to solve a system, but PARDISO is not as far as LAPACK can use global pivoting strategy. But it doesn't mean that LAPACK could solve everything - it can solve just a bit wider range of matrices, that's all.
I think it would be great if you send us your matrix in order we can experiment with it and try to help you.
Best regards,
Konstantin
0 Kudos
Danesh_Daroui
Beginner
2,118 Views
I expect my coefficient matrix to have very small eigenvalues which makes it very close to singular. Otherwise, it would be singular and LAPACK could also fail, that I know. There will be four vectors, three for the coefficient matrix in sparse format (non-zeros, column and row) and one for the right-hand-side. Would it be OK if I send these four vectors in MATLAB format so you can run same test on your side and see what will happen?

Thank you so much in advance,

D.
0 Kudos
Danesh_Daroui
Beginner
2,118 Views
Well, Now I am sure that PARDISO has a bug. I could successfully solve my problem using MUMPS package. What I can say is PARDISO has a problem (probably in pivoting) when the matrix is unsymmetric and complex and highly indefinite. I am sure about it. This problem (incorrect results) happened even with PARDISO solver from its original developers but MKL PARDISO is even worse and it even return different results each time it is called which proves that the routine is not thread safe. I have submitted a bug report.

D.
0 Kudos
Sergey_Solovev__Inte
New Contributor I
2,118 Views
Danesh,
OK, send these four vectors for us and we will try to reproduce the problem.
0 Kudos
Danesh_Daroui
Beginner
2,118 Views
Yes sure. I have written a small program that reads and solves a small equation where you can see the incorrect results. I hope this would help. I will send it to the support ticket that I had submitted before.

Thanks for your help.

D.
0 Kudos
Konstantin_A_Intel
2,118 Views
Hi Danesh,
I have analyzed your matrix passing it through LAPACK's zgesvd routine in order to compute singular values (fortunately, 266x266 size allows to do it easily).
The conclusion is that the matrix is not only ill-conditioned, but even singular. Please look at the statistics for this matrix:
Singular values: max=133.488855, min=0.000000
Condition number = inf
All condition numbers:
133.488855 73.293987 2.725539 2.679350 2.641821 2.625590 2.608541 2.557884 2.538276 2.515709.
2.488354 2.471750 2.425410 2.419850 2.403760 2.396062 2.387028 2.385062 2.368260 2.323211.
2.303360 2.296533 2.274010 2.241106 2.229485 2.216525 2.176490 2.161701 2.145840 2.115991.
2.108372 2.086315 2.051317 2.039938 2.032053 2.030238 1.995768 1.983408 1.968421 1.963458.
1.954383 1.942443 1.909278 1.899220 1.867873 1.850786 1.831537 1.820993 1.811035 1.802728.
1.798178 1.785756 1.774059 1.766122 1.741148 1.731737 1.710596 1.695952 1.640719 1.630052.
1.621706 1.617739 1.589568 1.574067 1.558955 1.524820 1.509524 1.461670 1.439147 1.404993.
1.387764 1.350588 1.336505 1.302357 1.300369 1.289759 1.281100 1.269441 1.239696 1.208933.
1.192244 1.185149 1.168707 1.131820 1.092711 1.076464 1.044913 1.033550 1.019720 1.000381.
0.974200 0.932945 0.928091 0.883230 0.861114 0.839156 0.818413 0.790118 0.742767 0.661282.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000.
That means that about 2/3 of your equations in the matrix are linearly dependent and the system of equations has infinite number of solutions.
PARDISO solver IS NOT intended for work with SINGULAR matrices.
My recommendation for you is to try to use "LAPACK-> Linear Least Square (LLS) Problems" component of MKL. It allows solving rank-deficient problems - precisely as represented by your matrix.
Best regards,
Konstantin
0 Kudos
Danesh_Daroui
Beginner
2,118 Views
Hi,

Thanks for your anwer. Maybe this is a special case in my problems, but in all cases, LAPACK can solve the problem successfully. Besides, MUMPS, the other sparse direct solver can always solve my problem without any problem. Again, if the coefficient matrix is ill-conditioned or not, it has absolutely nothing to do with PARDISO because it is a direct solver and should solve the equation anyway. Besides, Intel's implementation of PARDISO is not thread-safe. This means that I get different results when I run my problems. Sometimes I get "zero pivot" error and sometimes not. I am 100% sure that this is a bug in MKL PARDISO. I explained my reasons in the premier support.

Regards,

D.
0 Kudos
Konstantin_A_Intel
2,119 Views
Hi Danesh,

It seems it was an error in my computations in the previous post - I passed mistakenly a bit changed matrix to SVD routine. I'm sorry! In fact, there're no zero singular values in the matrix and hence it's NOT singular! Condition number is quite big (~10^14). However, it should be enough to produce more or less sensible solution.
We'll look closely into the problem and keep you informed when any results of investigation will be ready.
Thank you very much for this report!
Regards,
Konstantin
0 Kudos
Danesh_Daroui
Beginner
2,118 Views
Hi Konstantin,

Thanks for your further tests. I know that my matrices are ill-conditioned and I have even created good preconditioner for iterative solvers, but as you mentioned, direct solvers (e.g. PARDISO, MUMPS, etc.) should be able to solve the equation anyway. As a guess, I think the bug is in the pivoting, because even when I turn on the scaling vectors and weighted matching switches, I get zero-pivot although I should never get zero pivot when these switches are on. Besides, I get different results from time to time which means the routine is not thread-safe.

Regards,

Danesh
0 Kudos
Konstantin_A_Intel
2,118 Views
Hi Danesh,
I found a bug leading to incorrect solution of your matrix. All works fine with the fix.
I'll notify you when the fix will be officially available.
Best regards,
Konstantin
0 Kudos
miramin
Beginner
2,118 Views
Hi Danesh,
I found a bug leading to incorrect solution of your matrix. All works fine with the fix.
I'll notify you when the fix will be officially available.
Best regards,
Konstantin

I think I have the same problem as Danesh.

Could please let me know when the fix is available.

The sooner the better!!


Thank u so much.

AMin
0 Kudos
Gennady_F_Intel
Moderator
2,118 Views
The fix available in the latest versions 10.3 Update 5 and 6 available now from registration center.
Please check how it will work on your side and let us know.
--GEnnady
0 Kudos
Reply