- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
HI!
I used MKL's pardiso to get the decomposition of one of my LDLT matrices. I tried many methods, and even modified it directly in example pardiso_sym_getdiag.c, and finally reported an error in 332, unable to get the D matrix in LDL, the following is my The code, where the Ka matrix is the matrix obtained from the previous calculation, I have used mkl_sparse_s_export_csr to export the Ka matrix to export the triple and save it as a TXT file.
thank you for your help!
mkl_sparse_s_export_csr(Ka, &indexing, &rows, &cols, &row_strat, &rows_end, &col_index, &values);
a = values;
ia = row_strat;
ja = col_index;
MKL_INT n = rows;
MKL_INT mtype = -2; /* Real symmetric matrix */
/* RHS and solution vectors. */
float b[462], invdiag[462];
MKL_INT nrhs = 1; /* Number of right hand sides. */
/* Internal solver memory pointer pt, */
/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
/* or void *pt[64] should be OK on both architectures */
void* pt[64];
/* Pardiso control parameters. */
MKL_INT iparm[64];
MKL_INT maxfct, mnum, phase, error, msglvl;
/* Auxiliary variables. */
MKL_INT i;
float ddum; /* Double dummy */
MKL_INT idum; /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
for (i = 0; i < 64; i++)
{
iparm[i] = 0;
}
iparm[0] = 1; /* No solver default */
iparm[1] = 2; /* Fill-in reordering from METIS */
iparm[3] = 0; /* No iterative-direct algorithm */
iparm[4] = 0; /* No user fill-in reducing permutation */
iparm[5] = 0; /* Write solution into x */
iparm[6] = 0; /* Not in use */
iparm[7] = 0; /* Max numbers of iterative refinement steps */
iparm[8] = 0; /* Not in use */
iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
iparm[11] = 0; /* Not in use */
iparm[12] = 1; /* Maximum weighted matching algorithm is switched-on (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */
iparm[13] = 0; /* Output: Number of perturbed pivots */
iparm[14] = 0; /* Not in use */
iparm[15] = 0; /* Not in use */
iparm[16] = 0; /* Not in use */
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 */
iparm[34] = 1;
maxfct = 1; /* Maximum number of numerical factorizations. */
mnum = 1; /* Which factorization to use. */
msglvl = 0; /* Do not print statistical information in file */
error = 0; /* Initialize error flag */
/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
for (i = 0; i < 64; i++)
{
pt[i] = 0;
}
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
phase = 11;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
if (error != 0)
{
printf("\nERROR during symbolic factorization: " IFORMAT, error);
exit(1);
}
printf("\nReordering completed ... ");
printf("\nNumber of nonzeros in factors = " IFORMAT, iparm[17]);
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
phase = 22;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
if (error != 0)
{
printf("\nERROR during numerical factorization: " IFORMAT, error);
exit(2);
}
printf("\nFactorization completed ... ");
/* -------------------------------------------------------------------- */
/* .. Get inverse to the diagonal matrix D. Iterative refinement must be*/
/* turn off. */
/* -------------------------------------------------------------------- */
iparm[7] = 0;
phase = 332;
for (i = 0; i < n; i++)
{
b[i] = 1.0;
}
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, invdiag, &error);
if (error != 0)
{
printf("\nERROR during solution: " IFORMAT, error);
exit(3);
}
printf("\n Number of pertubed pivots = " IFORMAT " \n", iparm[13]);
if (iparm[13] != 0)
{
printf("\n The diagonal value were pertubed \n");
exit(4);
}
/* -------------------------------------------------------------------- */
/* .. Inverse elements of the array invdiag which are inversed diagonal */
/* values D^(-1) and print them . */
/* -------------------------------------------------------------------- */
printf("\n Print values of D for LDLT decomposition \n");
for (i = 0; i < n; i++)
{
invdiag[i] = 1.0 / invdiag[i];
printf("\n invdiag [" IFORMAT "] = %30.15e", i, invdiag[i]);
}
printf("\n");
/* -------------------------------------------------------------------- */
/* .. Termination and release of memory. */
/* -------------------------------------------------------------------- */
phase = -1; /* Release internal memory. */
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, &ddum, ia, ja, &idum, &nrhs,
iparm, &msglvl, &ddum, &ddum, &error);
- Tags:
- pardiso
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Huang,
Since your pardiso error value returns -1, it indicates that there is inconsistency in the input. You can get Pardiso to reveal more information about the inconsistencies by setting iparm[26] = 1 before calling Pardiso.
This setting helps you to see some more useful information.
For more information regarding iparm[26], please refer the below link:
Regards,
Vidya.
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Huang,
Thanks for reaching out to us.
I just checked with pardiso_sym_getdiag.c example code and it is working fine without any errors.
>> where the Ka matrix is the matrix obtained from the previous calculation
It would be a great help if you can share the complete reproducer code and the steps you have followed to run it so that we can check it from our end.
Please do let us know your OS details and MKL version as well.
Regards,
Vidya.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Huang,
Since your pardiso error value returns -1, it indicates that there is inconsistency in the input. You can get Pardiso to reveal more information about the inconsistencies by setting iparm[26] = 1 before calling Pardiso.
This setting helps you to see some more useful information.
For more information regarding iparm[26], please refer the below link:
Regards,
Vidya.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi,
Reminder:
Could you please let us know if your issue is resolved and provide us with an update regarding the same?
Regards,
Vidya.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you very much for your help!
I used double precision instead of single precision, and achieved good results, probably single precision is not suitable for my calculation problem.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Huang,
Thanks for the confirmation and glad to know that your issue is resolved.
As the issue is resolved we are closing this thread. Please post a new question if you need any additional information from Intel as this thread will no longer be monitored.
Have a Great Day!
Regards,
Vidya.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page