problem with multi-threading in pardiso

jjforest · ‎01-31-2009

I am not able to use multi-threading in pardiso. The following is the code I used in the program. In the solver progress message, I always get to use only thread 0. such as "In thread 0, at stage Pardiso: factor, steps passed (some percentage)". The matrix I used is huge and pardiso should use multi-threaded solver in my dual-core macbook pro (vista os). I am linking the program to pardiso using " mkl_intel_c.lib mkl_intel_thread.lib mkl_core.lib libguide40.lib". What else do I need to do in order to make use of multi-threaded pardiso solver.

Thanks,

jj

// testIntelParallel.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"

#include
#include
#include
#include
#include
#include
#include
using namespace std;

// libguide40.lib is for dynamically linked OpenMP (threading)

/* PARDISO prototype. */
#if defined(_WIN32) || defined(_WIN64)
#define pardiso_ PARDISO
#else
#define PARDISO pardiso_
#endif
#if defined(MKL_ILP64)
#define MKL_INT long long
#else
#define MKL_INT int
#endif

extern "C" int mkl_get_max_threads();
extern "C" void MKL_Set_Num_Threads( int number );

extern "C" MKL_INT PARDISO
(void *, MKL_INT *, MKL_INT *, MKL_INT *, MKL_INT *, MKL_INT *,
double *, MKL_INT *, MKL_INT *, MKL_INT *, MKL_INT *, MKL_INT *,
MKL_INT *, double *, double *, MKL_INT *);

#define BUFLEN 16
extern "C" int MKL_PROGRESS( int* ithr, int* step, char* stage, int len )
{

char buf[BUFLEN];

if( len >= BUFLEN ) len = BUFLEN-1;

strncpy( buf, stage, len );

buf[len] = '\0';

printf( "In thread %i, at stage %s, steps passed %i\n", *ithr, buf, *step );

return 0;

}

void test(int& nEqns, double*& values, int*& rowIndex, int*& columns, double*& rhs)
{
std::ifstream fileIn( "c:\\test.m", std::ios::binary );
if(std::ios::failbit & fileIn.rdstate()) {
return;
}

fileIn.read( (char *) &nEqns, sizeof(int) );
rowIndex = new int[nEqns + 1];
fileIn.read( (char *) rowIndex, sizeof(int) * (nEqns + 1));

int nNonzeros = rowIndex[nEqns];
values = new double[nNonzeros];
columns = new int[nNonzeros];
rhs = new double[nEqns];

fileIn.read( (char *) values, sizeof(double) * nNonzeros);
fileIn.read( (char *) columns, sizeof(int) * nNonzeros);
fileIn.read( (char *) rhs, sizeof(double) * nEqns);
}

MKL_INT main( void ) {

MKL_INT nEqns;
MKL_INT nNonZeros;
double* a = 0;
int* ia = 0;
int* ja = 0;
double* b = 0;
double* x = 0;

printf("construct stiffness matrix\n");

test(nEqns, a, ia, ja, b);

MKL_INT n = nEqns;

x = new double[nEqns];

printf("number of equations = %d\n", nEqns);

printf("number of CPUs = %d\n", mkl_get_max_threads());
MKL_Set_Num_Threads(mkl_get_max_threads()); // multiple threads

//MKL_INT mtype = -2; /* Real symmetric matrix */
MKL_INT mtype = 2; /* Positive Definite */
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 = 0;
}
iparm[0] = 1; /* No solver default */
iparm[1] = 2; /* Fill-in reordering from METIS */
/* Numbers of processors, value of MKL_NUM_THREADS */
iparm[2] = mkl_get_max_threads();
iparm[3] = 0; /* 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] = 2; /* Max numbers of iterative refinement steps */
iparm[8] = 0; /* Not in use */
iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
iparm[11] = 0; /* Not in use */
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[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] = -1; /* Output: Mflops for LU factorization */
iparm[19] = 0; /* Output: Numbers of CG Iterations */
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 = 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]);

printf("\n************************************************************************************");
DWORD t1 = timeGetTime();

/* -------------------------------------------------------------------- */
/* .. 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\n", error);
exit(2);
}
printf("\nFactorization completed ... ");

/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
phase = 33;
iparm[7] = 2; /* Max numbers of iterative refinement steps. */

FILE* fp;
fp=fopen("c:\\sol.txt","w");

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\n", error);
exit(3);
}

DWORD t2 = timeGetTime();
printf("\nSolve completed ... ");
printf("\n************************************************************************ time spent %d ... ", t2 - t1);
printf("\nThe solution of the system is: ");

for (i = 0; i < n; i++) {
fprintf(fp, "\n x [%d] = %g", i, x );
}
printf ("\n");

fclose(fp);

/* -------------------------------------------------------------------- */
/* .. 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);

delete []a;
delete []ia;
delete []ja;
delete []b;
delete []x;

return 0;
}

Konstantin_A_Intel · ‎07-16-2009

Hi jj,

We aresorry for the late reply.

I believe there're no problems with threading into your application. In fact, progress routine returns information only for the master thread, but it does not mean that PARDISO uses only a single thread inside.

In order to be sure, please set msglvl=1 into your main program. You should find a string like this one into your output:

< Parallel Direct Factorization with #processors: > 4

This is an actual number of threads used in PARDISO.

Regards,
Konstantin