topic it should be iparm[30]=1 on in Intel® oneAPI Math Kernel Library

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

Xiaofeng_W_ — Thu, 21 Jul 2016 02:24:33 GMT

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]);
  }
}

Hi,

Alexander_K_Intel2 — Thu, 21 Jul 2016 04:14:45 GMT

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

it should be iparm[30]=1 on

Xiaofeng_W_ — Thu, 21 Jul 2016 04:18:00 GMT

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

the issue is escalated and

Gennady_F_Intel — Sun, 24 Jul 2016 06:51:22 GMT

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

Hi,

Alexander_K_Intel2 — Fri, 29 Jul 2016 05:16:01 GMT

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