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 = 1; /* No solver default */
358 iparm = 2; /* Fill-in reordering from METIS */
359 iparm = 8;
360 iparm = 0; /* CGS */
361 iparm = 0; /* No user fill-in reducing permutation */
362 iparm = 0; /* Write solution into x */
363 iparm = 0; /* Not in use */
364 iparm = 0; /* Max numbers of iterative refinement steps */
365 iparm = 0; /* Not in use */
366 iparm = 13; /* Perturb the pivot elements with 1E-13 */
367 iparm = 1; /* Use nonsymmetric permutation and scaling MPS */
368 iparm = 0; /* Not in use */
369 iparm = 0; /* Not in use */
370 iparm = 0; /* Output: Number of perturbed pivots */
371 iparm = 0; /* Not in use */
372 iparm = 0; /* Not in use */
373 iparm = 0; /* Not in use */
374 iparm = -1; /* Output: Number of nonzeros in the factor LU */
375 iparm = -1; /* Output: Mflops for LU factorization */
376 iparm = 0; /* Output: Numbers of CG Iterations */
377 iparm = 1; /* check the data structure */
378 iparm = 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
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.
To have Pardiso check your matrix, set iparm=1 (not iparm=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
mecej4 wrote: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.
These make no sense. It is the indices of arrays b and x that may be zero- or one-based, not their values.
mecej4 wrote:This is done by m_row.push_back(m_val.size()+1);
the last entry should be set as m_row[m_dim]=m_nnz+1
mecej4 wrote:I have just tried this... My matrix seems to be correctly input.
To have Pardiso check your matrix, set iparm=1
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 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.