How to perform Schur Factorization and Forward/Backward properly with Pardiso？

[Pardiso_user]

Hi all, I'm new to oneMKL Pardiso. I want to factor a real unsym matrix $A$ as  ,

then I need to perform

forward substitution, i.e.,  to solve linear systems with triangular matrix $[L11,  0 ；L12， I ]$ 

and backward substitution, i.e., to solve linear systems with triangular matrix $[U11, U21; 0, I ]$.

I tried to set iparm[35] = -2 and iparm[23] = 1 to compute $S$, and phase = 331 to perform forward, phase = 333 to perform backward.

It seems that the forward substitution phase is good, but the backward substitution result is bad.

For backward substitution, I expect that the last components (associated with the schur part) of the solution are the same as the last components of the input, but they are not.

Could you please kindly provide an example ?

[Fengrui]

Did you perform phase=332? There is an example to calculate Schur complement and use it to perform the subsequent solve stages. The example is in the oneMKL package and located at $MKLROOT/share/doc/mkl/examples/examples_core_c.tgz. The source file is pardiso_schur.c in the sparse_directsolvers folder.

[Pardiso_user]

Hi Fengrui.  Thanks for your advice.

I learned the examples and I am still confused.

Here is my demo code.

I expect that the last two components of $x$ are the same as the last two components of $b$ and equal to one, but the results are not.

I get x[6] = { 0.45, 1.311688, 0.999091, -0.45, -0.909091, 0.11}.

I wonder if there is a mistake in my code.

Thanks !

// OS：Centos 7.7

MKL_INT n = 6;

MKL_INT ia[7] = {0, 2, 4, 7, 9, 11, 13};

MKL_INT ja[13] = {0, 5, 1, 4, 2, 3, 4, 3, 5, 4, 5, 3, 5};

double a[13] = {1, 5, 7, 9, 5, 1, 5, -1, 5, 11, 100, 1, 5};

MKL_INT mtype = 11;

MKL_INT nrhs = 1, maxfct = 1, mnum = 1, error = 0, msglvl = 0;

MKL_INT idum;

double ddum;

void *pt[64];

MKL_INT n_schur = 2;

int perm[6] = {0, 0, 0, 0, 1, 1};

pardisoinit(pt, &mtype, iparm);

iparm[34] = 1;

iparm[35] = -2;

iparm[23] = 1;

iparm[4] = 0;

iparm[30] = 0;

iparm[7] = 0;

phase = 11;

pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

int schur_nnz = iparm[35];

std::vector<int> iS(n_schur+1);

std::vector<int> jS(schur_nnz);

std::vector<double> S(schur_nnz);

int step = 1;

pardiso_export(pt, S.data(), iS.data(), jS.data(), &step, iparm, &error);

phase = 22;

iparm[35] = -2;

pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

phase = 333;

double x[6] = {0, 0, 0, 0, 0, 0};

double b[6] = {1, 1, 1, 1, 1, 1};

pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);