Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Pardiso: Schur complement for an SPD matrix

Zohar
Novice
378 Views

I have the following system: Ax = b, where A is a 2x2 block matrix:

A = [H   A1  
     A1' Z ]

H is SPD (symmetric positive definite) and Z is made of zeros.

(It's a Hessian with linear constraints.)

I'm solving the system using Schur complement (iparm[35]=2).

It's fine with mtype=-2 (symm, indefinite), but the factorization fails (-3, bad input) when mtype=2 (SPD).

Just to test the code, I tried it on an SPD A, which is solved fine with mtype=2.

It seems that mtype=2 requires that A is SPD, but for Schur it should only require H to be SPD?

0 Kudos
5 Replies
Zohar
Novice
319 Views

For a refresher on Schur complement, see section 3.5 in

https://igl.ethz.ch/projects/sparse-cholesky-update/sparse-cholesky-update-paper.pdf

My system is similar to eq. (20).

Only A is factorized, and it's SPD (mtype=2).

B is decomposed using lapacke (see the example pardiso_schur.c).

I can do the whole thing by myself: i) factorize A as SPD, ii) solve for B, iii) solve a small dense system with B. The bottleneck, however, is (ii). Forward/backward substitution for a single RHS is fast. For a matrix RHS, it's slow--which pardiso overcomes somehow.

0 Kudos
Shiquan_Su
Moderator
275 Views

Hi, Zohar:

Thanks for reporting these findings. Would you please provide more details of your test?

1, what is your OS?

2, what is your hardware?

3, what is your programming language?

4, would you please provide the source code of your test and the instruction of how you build your source code with Intel MKL?


0 Kudos
Zohar
Novice
262 Views

1. win11.

2. pc.

3. c++.

4. Use the mkl example "pardiso_schur.c". Feed to it the matrix:

A =
     1     0     1
     0     1     0
     1     0     0

where H = eye(2) and A1' = [1 0].
H is SPD and A is indefinte.

0 Kudos
Fengrui
Moderator
84 Views

Let M=[A B;C D], mtype is for the type of the whole matrix M. oneMKL Pardiso will decide internally how to compute for A^(-1)B when calculating the Schur complement of A which is D-C*A^(-1)*B. The user doesn't need to (actually can not) specify the type of A.

0 Kudos
Zohar
Novice
49 Views

Yeah, that was the point. M has a special structure where A is SPD while M isn't. Computing inv(A) using Cholesky (for SPD) is more efficient than L'DL (for symmetric only). Thus, it's a design question: shouldn't you let the user specify the type of A, and guessing the algorithm you use, it would be more efficient?

 

Specifically, if A is inverted, then why care about M definiteness at all?

0 Kudos
Reply