PARDISO. SCHUR Complement matrix

PARDISO. SCHUR Complement matrix

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

Zone: 

Thread Topic: 

Question
27 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

Hi,

Please see my comments below

Thanks,

Alex

- 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.

[akalinki] Current version of MKL Pardiso support schur complement matrix in dense format only

- 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.

[akalinki] You are correct, in current version of MKL you can handle this matrix via solving step but cannot access it's internal representation

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

[akalinki] file pardiso_schur_c in general mkl example folder show hot to calcluate schur complement of small matrix 

Thanks a lot.

Marc

I am sorry.  I am confused. I understand the Schur complement is a full matrix. This is fine, but I do not understand your answer about triangular matrices. You tell it can be used in the solving step. But what is the solving step doing when the Schur complement is used ? Solving for the full matrix A.X=B or something else. If the solving is done for the full matrix, it means that you need to factorize the Schur complement matrix in the solving step. Am I correct? Also, there is a remark in the documentation about using LAPACK to compute the Schur complement vector for large matrices. What do you call Schur complement vector? Is this the triangular matrices? Or the solution of the system or something else.

Sorry for asking may be naive questions.

Thanks

 

Quote:

marcsolal wrote:

I am sorry.  I am confused. I understand the Schur complement is a full matrix. This is fine, but I do not understand your answer about triangular matrices. You tell it can be used in the solving step. But what is the solving step doing when the Schur complement is used ? 

In therms from this paper phase 331 solve system with first, lower-triangular matrix, phase 333 with third, upper-triangular matrix

Quote:

marcsolal wrote:

 If the solving is done for the full matrix, it means that you need to factorize the Schur complement matrix in the solving step. Am I correct? 

Yes and no. If you set iparm[35] to 1 then only schur complement will be computed. If you set it to 2 then pardiso on factorization step will compute all factorized matrix and also will factorized Schur complement matrix. 

Quote:

marcsolal wrote:

Also, there is a remark in the documentation about using LAPACK to compute the Schur complement vector for large matrices. What do you call Schur complement vector? Is this the triangular matrices? Or the solution of the system or something else.

It is part of vector full size vector that correspond to Schut complement matrix. Actually that's the element of array x for which correspond positions in perm array set to 1.

Thanks for the answers. I think I need to be more specific. I have a (sparse) linear system with a 4 block matrix and I am trying to eliminate some dof from the system.

Let say my matrix is   A=[A11 A12;A21 A22]

My right hand side vector is B=[B1;B2]

And my system is:   A11 X1+A12X2=B1

                               A21 X1+A22X2=B2

I would like to express the system only in function of X2

I get: ( A22- A21A11-1A12)X2=B2-A21A11-1B1

The matrix on the left is the Schur complement matrix and PARDISO should give this. My problem is to be able to update the right hand side. A21A11-1B1 can be find from the triangular matrices. But if I do not have access to the triangular matrices, I need to find another way.

I hope I made myself clear.

Hi. Looks like there is a trick how it could be calculated without internal matrix data. I attach the photo with my calculations, and I really sorry for quality of this photo - i have only my telephone and paper from my daughter album for drawing :) 

Attachments: 

AttachmentSize
Downloadapplication/pdf Full page photo.pdf49.92 KB

Now I get it. It was exactly the meaning of my original question. I wanted to know if I can use the triangular matrix independently of the Schur complement matrix. I did not understand the meaning of forward and backward substitution. So the Schur complement matrix is not factorized in steps 331 and 333, which is exactly what I need. I am assuming step 332 in case of Schur complement matrix both factorizes and solves the Schur complement matrix which is a full matrix. It is why the documentation recommends to use Lapack routines to do that.

Thanks I think I can do what I need now.

 

Hi there,

I am trying to calculate large Schur complement for a symmetric sparse matrix using MKL PARDISO pardiso_schur_c.c program.

My parameters are:

MKL_INT n = 604034;

MKL_INT nonz = 2063839;

If I use the last 500 elements in perm vector = 1 (dimensions of S), then the S matrix is symmetric. If I increase the number of elements, say 1000, I get unsymmetric S.

What is the content of S? Does it contain the Schur complement, or the upper triangle is different from the lower triangle? Can anybody help me to understand this? Thanks in advance.

Konstantin

Quote:

Konstantin K. wrote:

Hi there,

I am trying to calculate large Schur complement for a symmetric sparse matrix using MKL PARDISO pardiso_schur_c.c program.

My parameters are:

MKL_INT n = 604034;

MKL_INT nonz = 2063839;

If I use the last 500 elements in perm vector = 1 (dimensions of S), then the S matrix is symmetric. If I increase the number of elements, say 1000, I get unsymmetric S.

What is the content of S? Does it contain the Schur complement, or the upper triangle is different from the lower triangle? Can anybody help me to understand this? Thanks in advance.

Konstantin

Hello Konstantin,

Please take a look at the related topic http://software.intel.com/en-us/articles/intel-mkl-support-to-new-functionality-schur-complement, I suppose it'll be helpful for understanding how Schur complement approach is implemented in MKL Pardiso.
Regarding your question, on the output Schur complement is returned as a general dense matrix.

Best regards,
Maria

 

Hi Maria,

 

Thank you very much for your response. I understand that the Schur matrix is a dense matrix with dimensions (n_schur,n_schur). I’ve been to the site you’ve suggested but my problem is quite different.

I have a large symmetric sparse matrix (upper triangle only). The task is to compute Schur complement on A22 , where A22 has different size. I have a check for symmetry of the Schur matrix at the end of the program. I’ve noticed that when the dimension of A22 is 500 x 500, the resulting Schur matrix is symmetric. If I increase the dimension, say 800 x 800, the Schur matrix I get is non-symmetrical. My understanding (and experience) is that if A is symmetric and can be partitioned as

|A11   A12|, A11 and A22 are invertible, then the Schur matrix, S = A22 - A21(A11) -1A12 is always

|A21   A22 |  symmetric.

My puzzle is why am I getting unsymmetrical result when the size is greater than 500? Is this a sign that something went wrong? Any thoughts? Thanks.

 

Kind Regards,

Konstantin

Hi Konstantin,

Can you please provide a reproducer, so I'll be able to investigate your case carefully?

Best regards,
Maria

Hi Maria,

 

Thank you very much for your response and willingness to help. Much appreciated. Let me give you some background of my story so you can understand the problem better.

My example comes from genetics. We have a matrix, say A that contains coefficients of relationships between animals in a given pedigree (probabilities of common genes coming from parents). This matrix is not easy to compute because it is very large (10M x10M). However, we have a number of animals that are genotyped (their genome is scanned) and we are interested only in those animals.

Therefore, we can partition A as

|A11  A12|

|A21  A22| , where A22 contains the animals of interest. We have an algorithm that can calculate A22 relatively inexpensive. However, we need the inverse of A22. Therefore the work to get the A22_inverse is doubled. Luckily, we have a fast algorithm to compute the inverse of A for millions of animals. However, we need only the portion of A_inverse that contains the genotyped animas, e.g. A22_inverse. That’s where the Schur complement from Pardiso comes into the picture. I have the entire A_inverse and I want to calculate A22 part only. Then the test A22 * A22_inverse = I, where I is identity matrix will tell me whether I am doing the right thing.  

My problem is that I get Schur = A22_inverse as an unsymmetrical matrix, which contradicts the theory.

There is another puzzle I found yesterday. I copied the lower triangle of the Schur matrix to the upper part and I got the required result (A22 * A22_inverse = I). This suggests to me that the program keeps the Schur result in the lower part of the matrix, which contradicts to the description of the matrix in the manual.

Anyway, attached please find the following files:

Test_schur.cpp – testing program, so you can reproduce the results.

Hol.ain – the matrix A_inverse in row compressed sparse format. The first part of the file contains IA and the second JA and A.

You can disable the use of Blitz array library. It is there because I have a version of the test program that uses Blitz arrays instead of C arrays.

Sorry for the long story, but I am very excited about using this program. It will save us lots of computing time. However, I need a proof that I am doing the right thing. Once again, thanks.

Kind regards,

Konstantin

Attachments: 

AttachmentSize
Downloadtext/x-c++src test_schur.cpp11.27 KB
Downloadapplication/rar Ainv.rar12.75 MB

Konstantin,

Thank you very much for the answer!
I was able to reproduce your problem and I can see that there is indeed problem with a lower triangular part of the matrix, so I'm working on it.
As a workaround I suggest you to use only upper triangular part of the Schur complement matrix in your computations.
Thanks for the issue.

Best regards,
Maria

Hi Maria,

I am glad that you’ve confirmed my findings. Just a comment: I believe that the problem is with the upper triangle, not the lower. In my other test program I copied the lower to the upper triangle and the test A22*A22_inv = I worked.

Kind regards,

Konstantin

Hello Konstantin,

You're right, of course. My apologies for the mix up.

Best regards,
Maria.
 

Hi Maria and Konstantin,

I'm also having issues with the Schur complement matrix returned from Pardiso not being symmetric when the original (sparse) matrix is symmetric. It seems that one of of the triangular parts of the returned matrix is correct. I'm only passing the one triangle of the original matrix to Pardiso.

Is it guaranteed that the triangular (seemingly correct) part of the Schur complement is indeed correct?

I'm using mkl 11.3 update 3. Is this issue fixed in Update 4?

Best and thanks,

Jens

 

Hi Jens,

You are right, currently only one triangular returned by Schur complement is correct and the issue still exist in last version of MKL - MKL2017u1

Thanks,

Alex

I am willing to use the Schur functionality of PARADISO on symmetric matrices. Can you please confirm which triangle is the correct one. Am I understanding correctly that the lower triangle is the one to use?

Thanks,

Marc 

Hi Marc,

You are right, lower triangle of the output matrix is correct and should be used.

Best regards,
Maria

 

Thanks. I am assuming Schur is not working with DSS handle. Right?

Thanks,

Marc

the problem is fixed in MKL 2017 u2. This update is released at Feb 22. Please check how this fix works on your side. 

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

 

 

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

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[i] = 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[i] = 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[i] = 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[i] = x[i];

    }

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

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

/* .. Solving Schur complement. */

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

    for(i = 0; i < n_schur; i++) ipiv[i] = 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[i]);

    }

    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

 

 

 

 

 

Hi Maria,

 

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

Regards,

 

Konstantin

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

 

 

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

Leave a Comment

Please sign in to add a comment. Not a member? Join today