Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Xiaofeng_W_
Beginner
30 Views

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

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

 

0 Kudos
4 Replies
30 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

 

Xiaofeng_W_
Beginner
30 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

Gennady_F_Intel
Moderator
30 Views

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

30 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

Reply