#include <stdio.h>
// #include <mkl.h>

/*
    Using MKL 2021.2.0
*/

int test(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;
}