- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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?
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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?

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page