topic Quote:marcsolal wrote: in Intel® oneAPI Math Kernel Library

PARDISO. SCHUR Complement matrix

marcsolal — Fri, 28 Oct 2016 17:41:59 GMT

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

Hi,

Alexander_K_Intel2 — Fri, 28 Oct 2016 18:35:19 GMT

Hi,

Please see my comments below

Thanks,

Alex

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

[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

marcsolal — Fri, 28 Oct 2016 18:53:48 GMT

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:

Alexander_K_Intel2 — Thu, 03 Nov 2016 04:15:22 GMT

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

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.

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

marcsolal — Thu, 03 Nov 2016 14:50:35 GMT

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=[A₁₁ A₁₂;A₂₁ A₂₂]

My right hand side vector is B=[B₁;B₂]

And my system is: A₁₁ X₁+A₁₂X₂=B₁

A₂₁ X₁+A₂₂X₂=B₂

I would like to express the system only in function of X₂

I get: ( A₂₂- A₂₁A₁₁^-1A₁₂)X₂=B₂-A₂₁A₁₁^-1B₁

_{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.}A₂₁A₁₁^-1B_{1 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

Alexander_K_Intel2 — Thu, 03 Nov 2016 17:55:25 GMT

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

Now I get it. It was exactly

marcsolal — Thu, 03 Nov 2016 18:41:26 GMT

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,

Konstantin_K_2 — Wed, 07 Dec 2016 06:20:50 GMT

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:

MariaZh — Wed, 07 Dec 2016 08:52:17 GMT

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,

Konstantin_K_2 — Wed, 07 Dec 2016 22:58:32 GMT

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 A²² , where A²² 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 A²² 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

|A¹¹ A¹²|, A¹¹ and A²²are invertible, then the Schur matrix, S = A²² - A²¹(A¹¹) ^-1A¹² is always

|A²¹ A²² | 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,

MariaZh — Thu, 08 Dec 2016 09:10:03 GMT

Hi Konstantin,

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

Best regards,
Maria

Hi Maria,

Konstantin_K_2 — Thu, 08 Dec 2016 23:18:21 GMT

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

Konstantin,

MariaZh — Fri, 09 Dec 2016 10:05:33 GMT

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,

Konstantin_K_2 — Sun, 11 Dec 2016 22:09:07 GMT

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,

MariaZh — Mon, 12 Dec 2016 09:05:46 GMT

Hello Konstantin,

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

Best regards,
Maria.

Hi Maria and Konstantin,

Jens_E_ — Tue, 13 Dec 2016 22:49:39 GMT

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,

Alexander_K_Intel2 — Fri, 16 Dec 2016 04:12:10 GMT

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

marcsolal — Thu, 12 Jan 2017 21:17:30 GMT

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,

MariaZh — Thu, 12 Jan 2017 21:22:58 GMT

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

marcsolal — Thu, 12 Jan 2017 21:41:15 GMT

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

Thanks,

Marc