topic Re: pardiso_schur.c with pardiso_export in Intel® oneAPI Math Kernel Library

pardiso_schur.c with pardiso_export

Qigeng — Mon, 13 Dec 2021 10:28:13 GMT

Hi,

I try to combine pardiso_schur.c with pardiso_export, and i want to get the Schur complement matrix and also compute all factorization arrays.

After read the develop reference and sample code of pardiso_export, i set iparms like this:

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] = -2; /* Use Schur complement */ iparm[24-1] = 1; 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 */

After pardiso_export, the schur results are:

schur_values: -7.64499 -7.72346 19.0363 schur_ia: 0 1 3 scur_ja: 0 0 1 pardiso_schur.c: -7.644985e+00 -7.723464e+00 -7.723464e+00 1.903631e+01

the schur_ia and schur_ja are not correct.

And after phase 22, i check the result of phase 331:

0 -7.4 1 1 1 0.588926 -2.14191 2.08939 pardiso_schur.c: 0 0.0931208 1 0.142857 0.2 -0.326816 -2.98962 2.08939

it doesn't match either.

Do I set the iparm wrong or other ipram should be set?

@Kirill_V_Intel

Thanks,

Re: pardiso_schur.c with pardiso_export

Kirill_V_Intel — Mon, 13 Dec 2021 19:07:30 GMT

Hi!

Could you please provide full code which you are running?

From what you described, the usage seems to be correct. And as a quick experiment, do you get correct Schur complement for iparm[36] = 2 (dense) or -1 (sparse)?

Best,
Kirill

Re: pardiso_schur.c with pardiso_export

Qigeng — Tue, 14 Dec 2021 03:11:17 GMT

I provide the patch of my code on pardiso_schur.c

I tried iparm[36] = -1 for sparse, the schur complement result is the same and result of phase 331 is not correct:

# result of schur complement: schur_values: -7.64499 -7.72346 19.0363 schur_ia: 0 1 3 schur_ja: 0 0 1 pardiso_schur.c: -7.644985e+00 -7.723464e+00 -7.723464e+00 1.903631e+01 # result of phase 331: 1 0.6 1.85714 1 1 0.588926 -2.14191 2.08939 pardiso_schur.c: 0 0.0931208 1 0.142857 0.2 -0.326816 -2.98962 2.08939

and for my understanding, the schur_ia should be: 0 2 3, and schur_ja should be: 0 1 1

Thanks

Re:pardiso_schur.c with pardiso_export

Gennady_F_Intel — Tue, 14 Dec 2021 03:57:17 GMT

Zhang,

it looks like a real problem and we will escalate it. The thread will be updated if needed.

Re: Re:pardiso_schur.c with pardiso_export

Qigeng — Tue, 14 Dec 2021 04:37:17 GMT

Get it!

Thanks.

Re: pardiso_schur.c with pardiso_export

Kirill_V_Intel — Tue, 14 Dec 2021 05:15:48 GMT

While it has no impact on the observed incorrectness, was it intentional that in your patch you additionally set iparm[36]=-1 (as well as iparm[36-1]=-1 for Schur complement option) which additionally adds VBSR feature?

Thanks,
Kirill

Re: pardiso_schur.c with pardiso_export

Qigeng — Tue, 14 Dec 2021 05:29:03 GMT

I'm not meanful to set iparm[36]=-1. I add two results, for the first one, I don't set iparm[36]=-1.

And I noticed your reply " do you get correct Schur complement for iparm[36] = 2 (dense) or -1 (sparse)?", so I think I can try it.

As a result, the patch has "iparm[36]=-1".

Thanks.

Re: pardiso_schur.c with pardiso_export

Kirill_V_Intel — Tue, 14 Dec 2021 06:14:20 GMT

Thanks for the clarification!

Based on the quick analysis, for symmetric matrix types we do return a lower-triangular CSR for Schur which explains what you see for Schur indexing arrays (they are correct as far as I understand, but just for the lower triangular instead of the upper-triangular which you expected).

As for why the result of phase 331 is incorrect, currently no updates.

Best,
Kirill

Re: pardiso_schur.c with pardiso_export

Qigeng — Tue, 14 Dec 2021 08:19:21 GMT

I see.

It's reasonable that the format is lower-triangular CSR.

BTW, I test pardiso and find both lower-triangular and upper-triangular csr format are accepeted.

Thanks.

Re: pardiso_schur.c with pardiso_export

Chao_Y_Intel — Fri, 06 Feb 2026 04:51:15 GMT

Hi, We fix this issue starting from oneMKL 2025.2 release. can you please check with new MKL releases: https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html Thanks for your report of this problem.

thanks,
Chao

Re: pardiso_schur.c with pardiso_export

rbunger — Thu, 23 Apr 2026 15:24:56 GMT

Hi,

I'm also interested in the same topic: The provided example "pardiso_schur.c", but with the factorization of the Schur matrix to solve a reduced system (plus optionally the sparse Schur matrix itself via "pardiso_export", but I think the factorization is good enough). In my case, the Schur matrix will be big so it has to be sparse.

Can somebody please provide such an (full) example? I couldn't download "0001_pardiso_schur.pardiso_export.txt".

Thanks,

Rainer

Re: pardiso_schur.c with pardiso_export

Fengrui — Sat, 02 May 2026 00:04:39 GMT

Hi Rainer,

The file, 0001_pardiso_schur.pardiso_export.txt, is attached. Please give it a try.

Best,

Fengrui

Re: pardiso_schur.c with pardiso_export

rbunger — Mon, 04 May 2026 08:27:00 GMT

Hi Fengrui,

I have used the fie attached to your message to patch the original "pardiso_schur.c" and it works, so I can build and run it.

Many thanks & regards,

Rainer

Re: pardiso_schur.c with pardiso_export

rbunger — Mon, 04 May 2026 10:09:48 GMT

Hmm, yes it compiles but the code is not correct:

the Schur matrix is printed like a dense Schur matrix (but it should be sparse)
there is only a 331 phase, but no 332 phase
the iparm array is possibly not correctly indexed (Fortran indexing instead of C indexing?)
in ILP64 mode the code crashes, possibly also due to the wrong (dense) printing loop

Re: pardiso_schur.c with pardiso_export

rbunger — Mon, 04 May 2026 10:12:35 GMT

... a 331 phase, but no 332 and 333 phase

Re: pardiso_schur.c with pardiso_export

rbunger — Mon, 04 May 2026 13:10:22 GMT

I think I have fixed the example somehow. But I honestly don't really know what I have done. I had another Fortran example and I have merged it into the above attached code. The results are correct (LP64). Remaining problem: It crashes for ILP64...

Re: pardiso_schur.c with pardiso_export

rbunger — Mon, 04 May 2026 13:26:00 GMT

Yet another problem is that there is no pardiso_export_64() function...

Re: pardiso_schur.c with pardiso_export

rbunger — Mon, 04 May 2026 13:37:43 GMT

Now it also works for ILP64, attached...