- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
I'm trying to use MKL pardiso (version MKL 2021.2) to obtain the inverse of Cholesky factor L^-1 by solving Lx = I (dentity) using phase 331 (here A = L * L' and A is well-conditioned).
But I found what I retrieve seems to be not matching up with the true L^-1 (despite ordering P). It's a bit confusing and I'm wondering if there is some parameter resulting in this phenomenon.
( The procedure was originally intended for a Lanczos algorithm where the largest eigenvalue of a matrix like L^-1 * X * L^-T is needed )
I've attached an example of 4 * 4 matrix with code reproducing this issue. Any help would be appreciated. Thank you in advance !
Best Regards.
int testpardiso(void) {
int n = 4;
/*
A =
0.331829000000000 0 0 -0.000157649000000
0 0.033182900000000 0 -0.022505800000000
0 0 0.331829000000000 0.022470700000000
-0.000157649000000 -0.022505800000000 0.022470700000000 0.580879000000000
*/
int Ap[5] = {0, 4, 7, 9, 10};
int Ai[10] = {0, 1, 2, 3, 1, 2, 3, 2, 3, 3};
double Ax[10] = {0.33182890474490068, 0, 0, -0.00015764871901342671, 0.33182890470071952,
0, -0.02250575799193312, 0.33182890465891607, 0.022470695607439466,
0.58087859335405345,};
double eye[16] = {1.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 1.0};
double aux[16] = {1.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 1.0};
double Linv[16] = {0.0};
int perm[4] = {0, 1, 2, 3};
int maxfct = 1; // Maximum number of factors
int mnum = 1; // The matrix used for the solution phase
int mtype = 2; // Real and symmetric positive definite
int msglvl = 0; // Print information
int idummy = 0; // Dummy variable for taking up space
int phs12 = 12;
int phs331 = 331;
int phsfree = -1;
int error = 0;
void *pdsWorker[64];
memset(pdsWorker, 0, sizeof(void *) * 64);
int params[64] = {
1, /* Non-default value */ 3, /* P Nested dissection */ 0, /* Reserved */
0, /* No CG */ 2, /* No user permutation */ 0, /* No overwriting */
0, /* Refinement report */ 0, /* Auto ItRef step */ 0, /* Reserved */
8, /* Perturb */ 0, /* Disable scaling */ 0, /* No transpose */
0, /* Disable matching */ 0, /* Report on pivots */ 0, /* Output */
0, /* Output */ 0, /* Output */-1, /* No report */
0, /* No report */ 0, /* Output */ 0, /* Pivoting */
0, /* nPosEigVals */ 0, /* nNegEigVals */ 0, /* Classic factorize */
0, 0, 1, /* Matrix checker */
0, 0, 0,
0, 0, 0,
0, 1, /* 0-based solve */ 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, /* No diagonal */ 0,
0, 0, 0,
0, 0, 0,
0
};
pardiso(pdsWorker, &maxfct, &mnum, &mtype, &phs12, &n,
Ax, Ap, Ai, perm, &idummy, params, &msglvl,
eye, aux, &error);
printf("Permutation: \n");
for (int i = 0; i < n; ++i) {
printf("%d, ", perm[i]);
}
printf("\n");
if (error) {
printf("Error code: %d \n", error);
goto cleanup;
}
pardiso(pdsWorker, &maxfct, &mnum, &mtype, &phs331, &n,
Ax, Ap, Ai, &idummy, &n, params, &msglvl,
eye, Linv, &error);
if (error) {
printf("Error code: %d \n", error);
goto cleanup;
}
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
printf("%10.6e, ", Linv[i + n * j]);
}
printf("\n");
}
cleanup:
pardiso(pdsWorker, &maxfct, &mnum, &mtype, &phsfree, &n,
NULL, NULL, NULL, &idummy, &idummy, params, &msglvl,
NULL, NULL, &error);
return error;
}
Testing environment:
MacBook Pro Apple Silicon 8GB RAM
Pardiso from MKL 2021.2 built using Rosetta
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
It seems to be an issue on user side. The A matrix in the reproducer (lines 61-64) does not match the A matrix actually used in the code (lines 6-11). In particular, the (2,2) element is wrong. It is 0.033182900000000, while in the code it is 0.33182890470071952.
If use the same A as in the code and then multiply that by Linv, we do get essentially the identity matrix as expected - the (2,2) element of that matrix is now 1 instead of 0.0976.
Regards,
Ruqiu
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi,
Thanks for reaching out to us.
As per the documentation, Ax[] Array should contain the non-zero elements but in the provided code we could see that the array contains 0 entries as well which might be the issue with the mismatch of the results.
So please refer to the below link for more details regarding the arguments of the pardiso call.
You can also look at the sample codes which are present in the installed directory of mkl for better understanding.
Regards,
Vidya.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
Thank you for your prompt reply. I believe that the answer is helpful but I still have a few concerns.
The zero elements in Ax[] array are numerical zeros due to cancellation instead of structural ones (there are consecutive solves where Ax[] is all non-zero). I have also tried perturbing the zeros by 1e-10 and the result seems unchanged.
I'm looking more closely into the documentation and will be grateful if there is any more comment available.
Best regards,
Gwz
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi,
>>what I retrieve seems to be not matching up with the true L^-1
Could you please share with us the expected results and the results that you are getting so that we could verify them from our end?
Regards,
Vidya.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
Thank you again for your help and for clarity I'm using the output format from MATLAB.
The full data matrix A is given by
A =
0.3318 0 0 -0.0002
0 0.0332 0 -0.0225
0 0 0.3318 0.0225
-0.0002 -0.0225 0.0225 0.5809
and the output L^-1 I extracted from Pardiso ( by Phase 331 ) is
Linv =
1.7360 0 0 0.0005
0.0000 1.7383 0 0.0673
-0.0000 -0.0046 1.7383 -0.0674
0 0 0 1.3121
From my understanding I should have A * transpose(Linv) * Linv = I but
>> A * Linv' * Linv
ans =
1.0000 0.0000 0.0000 0.0000
-0.0000 0.0976 0.0024 -0.0351
0.0000 -0.0000 1.0000 -0.0000
-0.0000 -0.0000 0.0000 1.0000
seems still a bit far from identity.
Best Regards,
Gwz
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi,
Thanks for your patience.
We are able to reproduce your issue. we are working on it and we will get back to you soon.
Regards,
Vidya.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
It seems to be an issue on user side. The A matrix in the reproducer (lines 61-64) does not match the A matrix actually used in the code (lines 6-11). In particular, the (2,2) element is wrong. It is 0.033182900000000, while in the code it is 0.33182890470071952.
If use the same A as in the code and then multiply that by Linv, we do get essentially the identity matrix as expected - the (2,2) element of that matrix is now 1 instead of 0.0976.
Regards,
Ruqiu
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thanks for your reply!
I have noticed the issue and there's no problem now.
Best regards,
Gwz

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