Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6998 Discussions

## Pardiso 11.3 does not work for iparm[30]=2 for partial solve

Beginner
554 Views

I recently tested the partial solve option with iparm[30]=1. I found setting only one element of perm to be 1 works, but when I set more than one element of perm of be 1 does not work. I also tried to link the same code with the pardiso from http://www.pardiso-project.org/, it works properly. Here the my test code of a small 4x4 matrix.

```#include <math.h>
#include <mkl.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

const MKL_INT size = 4;
MKL_INT n = size;
MKL_INT ia[5] = {1, 3, 5, 6, 7};
MKL_INT ja[6] = {1, 2, 2, 4, 3, 4};
double a[6] = {1, 5, 42, 6, 3, 4};
double rhs[size] = {0, 1, 0, 1};
MKL_INT perm[size] = {0, 1, 0, 1};
MKL_INT mtype = 2;  // positive definite
MKL_INT nrhs = 1;   /* Number of right hand sides. */
MKL_INT maxfct = 1; /* Maximum number of numerical factorizations. */
MKL_INT mnum = 1;   /* Which factorization to use. */
MKL_INT msglvl = 2; // verbose >= 1 ? verbose - 1
MKL_INT error = 0;  /* Initialize error flag */
MKL_INT iparm[64];
void *pt[64];

void Solve(const int iparm_30, double *x) {
iparm[30] = iparm_30;
MKL_INT phase = 11;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, &a[0], &ia[0], &ja[0],
&perm[0], &nrhs, &iparm[0], &msglvl, NULL, NULL, &error);
phase = 22;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, &a[0], &ia[0], &ja[0],
&perm[0], &nrhs, iparm, &msglvl, NULL, NULL, &error);
phase = 33;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, &a[0], &ia[0], &ja[0],
&perm[0], &nrhs, iparm, &msglvl, &rhs[0], &x[0], &error);
phase = -1;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, &a[0], &ia[0], &ja[0],
&perm[0], &nrhs, iparm, &msglvl, NULL, NULL, &error);

}

MKL_INT main(void) {
double x[4][size];
memset(&x[0][0], 0, sizeof(double) * 4 * size);

for (int i = 0; i < 64; i++)
pt = nullptr;

for (int i = 0; i < 64; i++)
iparm = 0;

iparm[0] = 1;  // No solver default
iparm[1] = 2;  // 0=minimum degree ordering, 2=Fill-in reordering from METIS
iparm[2] = 1;  // Numbers of processors, value of OMP_NUM_THREADS
iparm[3] = 0;  // 62; // No iterative-direct algorithm
iparm[4] = 0;  // No user fill-in reducing permutation
iparm[5] = 0;  // Write solution into x
iparm[6] = 0;  // Not in use
iparm[7] = 0;  // Max numbers of iterative refinement steps
iparm[8] = 0;  // Not in use
iparm[9] = 8;  // 13; // Perturb the pivot elements with 1E-13
iparm[10] = 0; // 1; // Use nonsymmetric permutation and scaling MPS
iparm[11] = 0; // Not in use
iparm[12] = 0; // matchings for highly indefinite symmetric matrices
iparm[13] = 0; // Output: Number of perturbed pivots
iparm[14] = 0; // Not in use
iparm[15] = 0; // Not in use
iparm[16] = 0; // Not in use
iparm[17] = 1; // Output: Number of nonzeros in the factor LU
iparm[18] = 0; // no Output: Mflops for LU factorization
iparm[19] = 0; // Output: Numbers of CG Iterations
iparm[20] = 1; // pivoting method

for (int iparm_30 = 0; iparm_30 <= 3; ++iparm_30) {
Solve(iparm_30, &x[iparm_30][0]);
}

printf("iparm30 %9d %10d %10d %10d\n", 0, 1, 2, 3);
for (int i = 0; i < size; ++i) {
printf("x[%d] = %10lf,  %10lf,  %10lf,  %10lf\n", i,  x[0], x[1], x[2], x[3]);
}
}
```

4 Replies
Employee
554 Views

Hi,

Could you confirm that you same output?

iparm30         0            1                  2          3
x[0] =   0.312500,    0.000000,    0.000000,    0.312500
x[1] =  -0.062500,   -0.062500,        -inf,       -0.062500
x[2] =   0.000000,   -0.062500,        -nan,      0.000000
x[3] =   0.343750,    0.343750,         inf,        0.343750

Look's like the problem in iparm[30]=2

Thanks,

Alex

Beginner
554 Views

it should be iparm[30]=1 on the tittle, I was not able to change it. My output is

iparm30     0            1               2               3

x[0] = 0.312500, 0.000000, 0.000000, 0.312500

x[1] = -0.062500, 0.000000, -0.062500, -0.062500

x[2] = 0.000000, -0.187500, 0.000000, 0.000000

x[3] = 0.343750, 0.531250, 0.343750, 0.343750

running under Linux Mint 18 64bit version (parallel_studio_xe_2016_composer_edition_for_cpp_update1 is the name of the installer package downloaded from intel). I linked the program with option -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5

Moderator
554 Views

the issue is escalated and you will be inform then the fix will be available into a official build of MKL.

Employee
554 Views

Hi,

We found the issue in code that provided incorrect result. Based on Intel MKL documentation:

The permutation vector perm must be present in all phases of Intel MKL PARDISO software. At the reordering step, the software overwrites the input vector perm by a permutation vector used by the software at the factorization and solver step.

The perm array changes after first call of pardiso so it's better to initialize it's before each new pardiso loop. If you reorginize code in such manner results will provide correctly.

Thanks,

Alex