Need some advice on the parameters of MKL's Pardiso

[Antoine__A_]

Hi,

I sent a post a couple of weeks ago, but I would like to specify my problem.

I want to solve a linear system of type Ax = b, for a common matrix A ( nonsymmetric, not positive-definite...), and b, x vectors.

Pardiso perfectly works in most of cases, but for my set of input data (see attached), it doesn't work: the result I obtain is totally aberrant. This is weird since A is inversible and very well conditionned (4.9819).

Here are my iparms :

 357         iparm[0] = 1; /* No solver default */
 358         iparm[1] = 2; /* Fill-in reordering from METIS */
 359         iparm[2] = 8;
 360         iparm[3] = 0; /* CGS */
 361         iparm[4] = 0; /* No user fill-in reducing permutation */
 362         iparm[5] = 0; /* Write solution into x */
 363         iparm[6] = 0; /* Not in use */
 364         iparm[7] = 0; /* Max numbers of iterative refinement steps */
 365         iparm[8] = 0; /* Not in use */
 366         iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
 367         iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
 368         iparm[11] = 0; /* Not in use */
 369         iparm[12] = 0; /* Not in use */
 370         iparm[13] = 0; /* Output: Number of perturbed pivots */
 371         iparm[14] = 0; /* Not in use */
 372         iparm[15] = 0; /* Not in use */
 373         iparm[16] = 0; /* Not in use */
 374         iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
 375         iparm[18] = -1; /* Output: Mflops for LU factorization */
 376         iparm[19] = 0; /* Output: Numbers of CG Iterations */
 377         iparm[27] = 1; /* check the data structure */
 378         iparm[31] = 1; /* iterative solver*/
 380         maxfct = 1; /* Maximum number of numerical factorizations. */
 381         mnum = 1; /* Which factorization to use. */

The Package ID of mkl is : l_mkl_p_10.0.011

Thanks for any help you can give me
Antoine

ps. A is stored in CSR format. Its file is composed of 4 lines : the size of datas (row column val), the ranks of the first element of the row, the columns, and the values. Tell me if you want another format.

[mecej4]

I notice from your data file that you use zero-based indexing. For this convention, you need to set iparm(35)=1 (in C++, iparm[34]=1) when calling Pardiso. Did you overlook this point? Secondly, it is not clear what you mean by "totally aberrant". To a casual reader, any results that contains numbers that are not very large or extremely small may appear reasonable. Therefore, you must state what criteria led you to say "totally aberrant".

[Antoine__A_]

This is what I mean by aberrant. The values of my matrix goes from -10²⁴ to 10⁴³. Actually, I should get a solution with all the values around 10⁴. I send you the x I calculated with Matlab, and the one I calculated with Pardiso. Yes, I overlooked iparm(35). This is not a problem of indexing.

[Alexander_K_Intel2]

Hi Antoine, Could you send example - pseudocode of using PARDISO for this matrix to check correctness of using it. With best regards, Alexander Kalinkin

[Antoine__A_]

Please find attached a file describing how I use Pardiso. Thanks for your time

[mecej4]

There are a number of errors in your .CPP program, mainly with regard to the CSR storage format, and one or two errors with regard to passing correct values in the array iparm[].

To have Pardiso check your matrix, set iparm[26]=1 (not iparm[27]=1 as you did).

The array ia must have one more entry: the last entry should be set as m_row[m_dim]=m_nnz+1; this is the CSR storage convention. You left this item undefined but, because of the preceding error, Pardiso did no checking of the matrix. You have loops such as for (int i = 0; i < m_b.size(); ++i) { m_b += 1; and for (int i = 0; i < m_x.size(); i++) { m_x -= 1; These make no sense. It is the indices of arrays b and x that may be zero- or one-based, not their values.

[Antoine__A_]

mecej4 wrote:

These make no sense. It is the indices of arrays b and x that may be zero- or one-based, not their values.

You are right, this really make no sens. Actually, I haven't done that in my real program. I needed to change the program a little and I added that, I don't know why.

mecej4 wrote:

the last entry should be set as m_row[m_dim]=m_nnz+1

This is done by m_row.push_back(m_val.size()+1);

mecej4 wrote:

To have Pardiso check your matrix, set iparm[26]=1

I have just tried this... My matrix seems to be correctly input.

[Antoine__A_]

I solved my problem thanks to a user from the Pardiso's mailing list. That comes from iparm[12]. I don't know why yet, but I needed to switch on this parameter. In my manual of MKL, the only description I have is : "iparm(12) : This parameter is reserved for future use. Its value must be set to 0." ;) This must be an ancient version. Thanks for your time, Antoine

[Alexander_K_Intel2]

Antoine, I see in your example iparm[11]=iparm[12]=0, correspondingly transpose solving and matching. How did you change it to resolve your issue? With best regards, Alexander Kalinkin

[Antoine__A_]

It works with iparm[12]=1.

[Antoine__A_]

Sorry, it's iparm[13] since I'm in C indexing. And so I understand why it is so important

[mecej4]

I think that your choosing iparm[0] = 1 was the root cause of many of the problems.
Had you chosen the solver defaults for your matrix type, you would have had a working program quickly, which you could then tweak by modifying the iparm[] values to your satisfaction.
However, you chose to set iparm[0] to a non-zero value at the outset, which then shifted the responsibility to you for setting all 64 values correctly, taking into account the documentation's 1-based indexing and your C++ program's 0-based indexing.

[Alexander_K_Intel2]

So iparm[12] or iparm[13] ? :) Setting iparm[12] turn on matching that strongly recommend for non-symmetric case With best regards, Alexander Kalinkin