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
- Problem with PARDISO

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

Danesh_Daroui

Beginner

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

03-30-2011
08:29 AM

451 Views

Problem with PARDISO

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

Mx=b

and then convert it to CSR format as below:

int nnz=0;

for (int i=0; i

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

1 Solution

Konstantin_A_Intel

Employee

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

05-26-2011
04:55 AM

451 Views

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

Link Copied

19 Replies

Konstantin_A_Intel

Employee

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

03-30-2011
08:22 PM

451 Views

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

Konstantin_A_Intel

Employee

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

03-30-2011
09:54 PM

451 Views

As you mentioned you're to use LAPACK, you can call SVD for your matrix and compute conditional number:

Regards,

Konstantin

Danesh_Daroui

Beginner

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

03-31-2011
02:27 AM

451 Views

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.

Konstantin_A_Intel

Employee

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

03-31-2011
04:24 AM

451 Views

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

Danesh_Daroui

Beginner

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

03-31-2011
04:50 AM

451 Views

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.

Sergey_K_Intel1

Employee

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

03-31-2011
07:35 AM

451 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

Danesh_Daroui

Beginner

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

03-31-2011
08:05 AM

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

Konstantin_A_Intel

Employee

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

03-31-2011
09:39 AM

451 Views

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

Danesh_Daroui

Beginner

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

03-31-2011
10:21 AM

451 Views

Thank you so much in advance,

D.

Danesh_Daroui

Beginner

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

04-18-2011
10:50 AM

451 Views

D.

Sergey_Solovev__Inte

New Contributor I

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

04-20-2011
01:10 AM

451 Views

Danesh,

OK, send these four vectors for us and we will try to reproduce the problem.

OK, send these four vectors for us and we will try to reproduce the problem.

Danesh_Daroui

Beginner

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

04-20-2011
03:45 AM

451 Views

Thanks for your help.

D.

Konstantin_A_Intel

Employee

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

05-25-2011
03:10 AM

451 Views

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

Danesh_Daroui

Beginner

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

05-26-2011
01:43 AM

451 Views

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.

Konstantin_A_Intel

Employee

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

05-26-2011
04:55 AM

452 Views

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

Danesh_Daroui

Beginner

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

05-29-2011
08:36 AM

451 Views

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

Konstantin_A_Intel

Employee

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

06-05-2011
10:48 PM

451 Views

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

miramin

Beginner

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

09-13-2011
03:12 PM

451 Views

Quoting Konstantin Arturov (Intel)

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

Could please let me know when the fix is available.

The sooner the better!!

Thank u so much.

AMin

Gennady_F_Intel

Moderator

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

09-13-2011
07:50 PM

451 Views

Please check how it will work on your side and let us know.

--GEnnady

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.