Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6998 Discussions

## Retrieve inverse of the Cholesky factor from MKL pardiso

Beginner
1,429 Views

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

1 Solution
Moderator
1,162 Views

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

7 Replies
Moderator
1,390 Views

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.

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso.html

You can also look at the sample codes which are present in the installed directory of mkl for better understanding.

Regards,

Vidya.

Beginner
1,375 Views

Hello,

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

Moderator
1,354 Views

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.

Beginner
1,346 Views

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

Moderator
1,305 Views

Hi,

We are able to reproduce your issue. we are working on it and we will get back to you soon.

Regards,

Vidya.

Moderator
1,163 Views

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

Beginner
1,132 Views