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

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

*!=(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.

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

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

Link Copied

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

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

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

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.

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

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

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.

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

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

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

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.

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

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

Thank you so much in advance,

D.

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

D.

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

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

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

Thanks for your help.

D.

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

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

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.

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

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

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

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

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

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

Could please let me know when the fix is available.

The sooner the better!!

Thank u so much.

AMin

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

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