PARDISO. SCHUR Complement matrix

[marcsolal]

I would like to use PARDISO to compute the SCHUR complement. I am trying to understand the documentation and I have questions:

- I understand that the Schur complement matrix is obtained in the solution vector. I had a look on PARDISO 5.0 (not the Intel software) documentation and the SCHUR complement is returned as a sparse matrix. Is MKL PARDISO returning a full matrix or a sparse matrix ? Also which element is where in the matrix.

- I will probably need to access the full factorization ie the Schur complement but also the other matrices (L11, L21, U11,U12 from the documentation). I understand it is possible to compute these matrices, but how can we access those and in which form.

- do you have a code example of PARDISO for computation of the SCHUR complement and access to all the matrices.

Thanks a lot.

Marc

[Konstantin_K_2]

Hi Maria,

It's Konstantin again. I am still having problems with the Schur complement. If the matrix size is below 50K it works well, e.g

A22 * A22_inv = I (identity matrix). However, when the size is greater than 50K the product of the Schur times the original matrix does not give identity matrix. I am using the lower triangle of the Schur matrix. Any thoughts? 

Regards,

Konstantin

[MariaZh]

Hi Konstantin,

Could you please clarify how you compute A22_inv?
Also, it would be great if you provide a reproducer for this case.

Best regards,
Maria

[Konstantin_K_2]

Hi Maria,

I have again problems with schur complement. Below is the code I used.

With the second example, n = 10, nonz = 28, the Schur complement is calculated correctly.

With the first example, n = 9, nonz = 18. I get some weired matrix.

The Schur should be:

   1.83    0.50   0.00  -1.00  -0.67    0.00
  0.50    1.50   0.00  -1.00   0.00    0.00
  0.00    0.00   1.00   0.00   0.00    0.00
 -1.00   -1.00   0.00   2.00   0.00    0.00
 -0.67    0.00   0.00   0.00   1.33    0.00
  0.00    0.00   0.00   0.00   0.00    1.00

The result from the program is :

Schur matrix: 
0.00  1.83  0.50 -1.00 -0.66 0.00
0.00  0.50  1.50 -1.00  0.00 0.00
1.00  0.00  0.00  0.00  0.00 0.00
0.00 -1.00 -1.00  0.00  2.00 0.00
0.00 -0.66  0.00  0.00  0.00 1.33
0.00  0.00  0.00  1.00  0.00 0.00

Can you have a look at the code. I am doing something wrong, I guess. Thanks in advance.

The code:

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include "mkl.h"

MKL_INT main (void)

{

    /* Matrix data. */

  //      MKL_INT n = 8;

    MKL_INT n = 9;

    MKL_INT ia[10] = {1, 4, 7, 9, 13, 15, 16, 17, 18,19};

    MKL_INT ja[18]= {1, 3, 9, 2, 4, 8, 3, 9, 4, 5, 7, 8, 5, 7, 6, 7, 8, 9};

    double a[18] = {  1.50,   0.50,  -1.00,   1.50,   0.50,  -1.00,   1.50,  -1.00,   2.00,   0.50,  -1.00,  -1.00,   1.50,  -1.00,   1.00,   2.00,   2.00,   2.00};

  /*

    MKL_INT n = 10;

    MKL_INT ia[11] = {1, 8, 10, 13, 18, 22, 24, 26, 27, 28, 29};

    MKL_INT ja[28]= {1, 2,    4,  5,  6,  7,      9,

      2,        5,

         3, 4,      6,

            4,  5,  6,  7,          10,

                5,      7,   8,     10,

                    6,          9,

                        7,   8,

                             8,

                                9,

                                    10};

   double a[28] =  {2.5,  0.5,        0.5, -1.0,  0.5, -1.0,       -1.0,

         1.5,             -1.0,

  1.5,  0.5,       -1.0,

                     2.5,  0.5, -1.0, -1.0,             -1.0,

                           3.0,        0.5, -1.0,       -1.0,

                                 2.5,             -1.0,

                                       2.5, -1.0,

                                             2.0,

                                                   2.0,

                                                         2.0};

  */

    int matrix_order=LAPACK_ROW_MAJOR;

    char uplo = 'U';

    MKL_INT mtype = -2;   /* Real symmetric matrix */

    /* RHS and solution vectors. */

    //        double b[8], x[8];

    double b[10], x[10];

    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, info;

    /* Auxiliary variables. */

    MKL_INT i, j;

    double ddum;          /* Double dummy */

    MKL_INT idum;         /* Integer dummy. */

    /* Schur data */

    /*

     double schur[4] = {0.0, 0.0,

                        0.0, 0.0 };

     MKL_INT perm[8] = {0, 0, 0, 0, 0, 0, 1, 1};

    MKL_INT ipiv[2];

    */

    //  MKL_INT n_schur = 2; /* Schur complement solution size */

      double schur[36] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0,

                          0.0, 0.0, 0.0, 0.0, 0.0, 0.0,

                          0.0, 0.0, 0.0, 0.0, 0.0, 0.0,

                          0.0, 0.0, 0.0, 0.0, 0.0, 0.0,

                          0.0, 0.0, 0.0, 0.0, 0.0, 0.0,

  0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

   MKL_INT perm[9] = {0, 0, 0, 1, 1, 1, 1, 1, 1};

     MKL_INT ipiv[6];

    MKL_INT n_schur = 6; /* Schur complement solution size */

/* -------------------------------------------------------------------- */

/* .. Setup Pardiso control parameters. */

/* -------------------------------------------------------------------- */

    for ( i = 0; i < 64; i++ )

    {

        iparm = 0;

    }

    iparm[1-1] = 1;         /* No solver default */

    iparm[2-1] = 2;         /* Fill-in reordering from METIS */

    iparm[10-1] = 8;        /* Perturb the pivot elements with 1E-13 */

    iparm[11-1] = 0;        /* Use nonsymmetric permutation and scaling MPS */

    iparm[13-1] = 0;        /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */

    iparm[14-1] = 0;        /* Output: Number of perturbed pivots */

    iparm[18-1] = -1;       /* Output: Number of nonzeros in the factor LU */

    iparm[19-1] = -1;       /* Output: Mflops for LU factorization */

    iparm[36-1] = 1;        /* Use Schur complement */

    maxfct = 1;           /* Maximum number of numerical factorizations. */

    mnum = 1;             /* Which factorization to use. */

    msglvl = 1;           /* 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 = 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, perm, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

    if ( error != 0 )

    {

        printf ("\nERROR during symbolic factorization: %d", error);

        exit (1);

    }

    printf ("\nReordering completed ... ");

/* -------------------------------------------------------------------- */

/* .. Numerical factorization. */

/* -------------------------------------------------------------------- */

    phase = 22;

    pardiso(pt, &maxfct, &mnum, &mtype, &phase,

            &n, &a, &ia, &ja, perm, &nrhs,

            iparm, &msglvl, &ddum, &schur, &error);

    if ( error != 0 )

    {

        printf ("\nERROR during numerical factorization: %d", error);

        exit (2);

    }

    printf ("\nFactorization completed ... ");

    printf("\nSchur matrix: \n"); fflush(0);

    for(i=0; i<n_schur; i++)

    {

        for(j=0; j<n_schur; j++)

        {

            printf("%f ",schur[i*n_schur+j]); fflush(0);

        }

        printf("\n"); fflush(0);

    }

/* -------------------------------------------------------------------- */

/* .. Reduce solving phase. */

/* -------------------------------------------------------------------- */

    phase = 331;

    /* Set right hand side to one. */

    for ( i = 0; i < n; i++ )

    {

        b = 1.0;

    }

    pardiso(pt, &maxfct, &mnum, &mtype, &phase,

            &n, a, ia, ja, perm, &nrhs,

            iparm, &msglvl, b, x, &error);

    if ( error != 0 )

    {

        printf ("\nERROR during solution phase 331: %d", error);

        exit (331);

    }

    for ( i = 0; i < n; i++ )

    {

        b = x;

    }

    printf ("\nSolve phase 331 completed ... ");

/* -------------------------------------------------------------------- */

/* .. Solving Schur complement. */

/* -------------------------------------------------------------------- */

    for(i = 0; i < n_schur; i++) ipiv = 0;

    info = LAPACKE_dsytrf(matrix_order,uplo,n_schur,&schur,n_schur,ipiv);

    if(info != 0)

    {

        printf("info after LAPACKE_dsytrf = %d \n", info);

        return 1;

    }

    info = LAPACKE_dsytrs(matrix_order,uplo,n_schur,nrhs,&schur,n_schur,&ipiv,&b[n-n_schur],nrhs);

    if(info != 0)

    {

        printf("info after LAPACKE_dsytrs = %d \n", info);

        return 1;

    }

/* -------------------------------------------------------------------- */

/* .. Expansion solving phase. */

/* -------------------------------------------------------------------- */

    phase = 333;

    /* Set right hand side to x one. */

    pardiso(pt, &maxfct, &mnum, &mtype, &phase,

            &n, a, ia, ja, perm, &nrhs, 

            iparm, &msglvl, b, x, &error);

    if ( error != 0 )

    {

        printf ("\nERROR during solution: %d", error);

        exit (333);

    }

    printf ("\nSolve phase 333 completed ... ");

    printf ("\nSolve completed ... ");

    printf ("\nThe solution of the system is: ");

    for ( i = 0; i < n; i++ )

    {

        printf ("\n x [%d] = % f", i, x);

    }

    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);

    return 0;

}

Kind regards,

Konstantin

[Konstantin_K_2]

Hi Maria,

Did you have a chance to look at my schur problem?

Regards,

Konstantin

[Gennady_F_Intel]

with the version of Intel MKL 2018 u3 as well as with the newest 2019 beta, I see  this Schur complement:

Factorization completed ... 

Schur matrix: 

1.833333 0.500000 0.000000 -1.000000 -0.666667 0.000000 

0.500000 1.500000 0.000000 -1.000000 0.000000 0.000000 

0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 

-1.000000 -1.000000 0.000000 2.000000 0.000000 0.000000 

-0.666667 0.000000 0.000000 0.000000 1.333333 0.000000 

0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 

MKL_VERBOSE Intel(R) MKL 2018.0 Update 3 Product build 20180406 for Intel(R) 64 architecture Intel(R) Advanced Vector Extensions 2 (Intel(R) AVX2) enabled processors, Win 2.60GHz cdecl intel_thread

[MariaZh]

Hi Konstantin,

I'm sorry for the delayed reply!
I've checked you program and the output is correct for the small (N = 9) case with the latest MKL version. 
Can you please double check on your side and let us know what you find out? 

Best regards,
Maria