Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
The Intel sign-in experience is changing in February to support enhanced security controls. If you sign in, click here for more information.
6680 Discussions

Retrieve inverse of the Cholesky factor from MKL pardiso

Gwz
Beginner
637 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

0 Kudos
1 Solution
Ruqiu_C_Intel
Employee
370 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


View solution in original post

7 Replies
VidyalathaB_Intel
Moderator
598 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/spars...

 

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

 

Regards,

Vidya.

 

Gwz
Beginner
583 Views

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

VidyalathaB_Intel
Moderator
562 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.


Gwz
Beginner
554 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

 

 

VidyalathaB_Intel
Moderator
513 Views

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.


Ruqiu_C_Intel
Employee
371 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


Gwz
Beginner
340 Views

Thanks for your reply!

I have noticed the issue and there's no problem now.

 

Best regards,

Gwz

Reply