Решение разреженного матричного уравнения с помощью Пардизо, созданного с помощью броненосца

Я хотел использовать Pardiso-Solver из MKL вместо встроенного решателя разреженных матриц. Для этого я сначала создаю матрицу в броненосце, а затем преобразую ее в CSC-матрицу, которая помещается в решатель Pardiso. Мой код

int main()
{
/* Matrix data. */
MKL_INT n = 8;
arma::sp_mat single_mat = arma::sp_mat(n, n);
single_mat.diag(-2).fill(1);
single_mat.diag(-1).fill(-8);
single_mat.diag(0).fill(0);
single_mat.diag(1).fill(8);
single_mat.diag(2).fill(-1);
single_mat = single_mat / 0.0012;
std::cout << "Number of nonzeros: " << single_mat.n_nonzero << '\n';
std::cout << "Number of n_cols: " << single_mat.n_cols << '\n';
MKL_INT *ia = new MKL_INT[single_mat.n_cols+1];
MKL_INT *ja = new MKL_INT[single_mat.n_nonzero];
double *a = new double[single_mat.n_nonzero];
for (size_t i = 0; i < single_mat.n_nonzero; i++)
ja[i] = single_mat.row_indices[i];
for (size_t i = 0; i < single_mat.n_nonzero; i++)
a[i] = single_mat.values[i];
for (size_t i = 0; i < single_mat.n_cols+1; i++)
ia[i] = single_mat.col_ptrs[i];std::cout << "\n" << arma::mat(single_mat) << "\n";

MKL_INT mtype = 11;       /* Real symmetric matrix */
/* RHS and solution vectors. */
double *b = new double[n], *x = new double[n], *bs = new double[n], res, res0;
bs[0] = 0;
bs[1] = 0;
MKL_INT nrhs = 1;     /* Number of right hand sides. */
/* Internal solver memory pointer pt, */
/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
/* or void *pt[64] should be OK on both architectures */
void *pt[64];
/* Pardiso control parameters. */
MKL_INT iparm[64];
MKL_INT maxfct, mnum, phase, error, msglvl;
/* Auxiliary variables. */
MKL_INT i;
double ddum;          /* Double dummy */
MKL_INT idum;         /* Integer dummy. */
/* -------------------------------------*/
/* .. Setup Pardiso control parameters. */
/* -------------------------------------*/
for (i = 0; i < 64; i++)
{
iparm[i] = 0;
}
iparm[0] = 1;         /* No solver default */
iparm[1] = 2;         /* Fill-in reordering from METIS */
iparm[3] = 0;         /* No iterative-direct algorithm */
iparm[4] = 0;         /* No user fill-in reducing permutation */
iparm[5] = 0;         /* Write solution into x */
iparm[7] = 2;         /* Max numbers of iterative refinement steps */
iparm[9] = 13;        /* Perturb the pivot elements with 1E-13 */
iparm[10] = 1;        /* Use nonsymmetric permutation and scaling MPS */
iparm[12] = 0;        /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */
iparm[13] = 0;        /* Output: Number of perturbed pivots */
iparm[17] = -1;       /* Output: Number of nonzeros in the factor LU */
iparm[18] = -1;       /* Output: Mflops for LU factorization */
iparm[19] = 0;        /* Output: Numbers of CG Iterations */
iparm[26] = 1;        /* Checks if matrix is correct*/
iparm[34] = 1;        /* PARDISO use C-style indexing for ia and ja arrays */
maxfct = 1;           /* Maximum number of numerical factorizations. */
mnum = 1;         /* Which factorization to use. */
msglvl = 0;           /* Print statistical information in file */
error = 0;            /* Initialize error flag */
/* ----------------------------------------------------------------*/
/* .. Initialize the internal solver memory pointer. This is only  */
/*   necessary for the FIRST call of the PARDISO solver.           */
/* ----------------------------------------------------------------*/
for (i = 0; i < 64; i++)
{
pt[i] = 0;
}
/* --------------------------------------------------------------------*/
/* .. Reordering and Symbolic Factorization. This step also allocates  */
/*    all memory that is necessary for the factorization.              */
/* --------------------------------------------------------------------*/
phase = 11;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
if (error != 0)
{
printf("\nERROR during symbolic factorization: %d", error);
exit(1);
}
printf("\nReordering completed ... ");
printf("\nNumber of nonzeros in factors = %d", iparm[17]);
printf("\nNumber of factorization MFLOPS = %d", iparm[18]);
/* ----------------------------*/
/* .. Numerical factorization. */
/* ----------------------------*/
phase = 22;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
if (error != 0)
{
printf("\nERROR during numerical factorization: %d", error);
exit(2);
}
printf("\nFactorization completed ... ");
/* -----------------------------------------------*/
/* .. Back substitution and iterative refinement. */
/* -----------------------------------------------*/
phase = 33;
iparm[7] = 2;         /* Max numbers of iterative refinement steps. */
/* Set right hand side to one. */
for (i = 0; i < n; i++)
{
b[i] = 1;
}
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);
if (error != 0)
{
printf("\nERROR during solution: %d", error);
exit(3);
}
printf("\nSolve completed ... ");
printf("\nThe solution of the system is: ");
for (i = 0; i < n; i++)
{
printf("\n x [%d] = % f", i, x[i]);
}
printf("\n");
//Check with Lapack solve
arma::colvec b_vec = arma::colvec(n);
b_vec.fill(1);
arma::colvec x_vec = arma::zeros(n);
arma::spsolve(x_vec, single_mat, b_vec, "lapack");//End Lapack
printf("\nx_vec is: ");
for (i = 0; i < n; i++)
{
printf("\n x [%d] = % f", i, x_vec[i]);
}
printf("\n");
//Check resolutions
arma::colvec x_pardiso = arma::colvec(n);
for (size_t i = 0; i < n; i++)
x_pardiso[i] = x[i];

std::cout << "\n" << single_mat*x_pardiso << '\n' << single_mat*x_vec << '\n';

/* --------------------------------------*/
/* .. Termination and release of memory. */
/* --------------------------------------*/
phase = -1;           /* Release internal memory. */
PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
&n, &ddum, ia, ja, &idum, &nrhs,
iparm, &msglvl, &ddum, &ddum, &error);
std::cout << "\nPardiso freed memory\n";
delete[]a;
delete[]ja;
delete[]ia;
delete[]x;
delete[]b;
delete[]bs;
return 0;
}

Я проверяю результат, используя Lapack-сольвер из броненосца. Но в то время как последний дает мне правильный результат, я не получаю правильный результат от pardiso-solver. Мой текущий вывод

The solution of the system is:
x [0] =  0.000674
x [1] = -0.000089
x [2] =  0.000488
x [3] = -0.000287
x [4] =  0.000288
x [5] = -0.000487
x [6] =  0.000089
x [7] = -0.000677

x_vec is:
x [0] = -0.000673
x [1] =  0.000089
x [2] = -0.000487
x [3] =  0.000287
x [4] = -0.000287
x [5] =  0.000487
x [6] = -0.000089
x [7] =  0.000673

Значения для x_vec являются правильными, значения для x — нет. Зачем?