pardiso_schur.c with pardiso_export

[Qigeng]

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,

[Kirill_V_Intel]

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

[Qigeng]

Hi

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

[Gennady_F_Intel]

Zhang,

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

[Qigeng]

Get it!

Thanks.

[Kirill_V_Intel]

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

[Qigeng]

Hi

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.

[Kirill_V_Intel]

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

[Qigeng]

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.

[Chao_Y_Intel]

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

[rbunger]

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

[Fengrui]

Hi Rainer,

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

Best,

Fengrui

[rbunger]

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

[rbunger]

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

[rbunger]

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

[rbunger]

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

[rbunger]

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

[rbunger]

Now it also works for ILP64, attached...