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.