Pardiso : memory release problem (MKL 10.1.1.019)

[minos67paran_com]

Hi, I am trying to use the Pardiso solver for the nonstationary Navier-Stokes equations.
The only problem is that the system memory is increasing constantly.

I tried MKL_FreeBuffers(), but it was not working.
I am curious to know if the problem could be fixed.
Is anyone who use Pardiso solver to the dynamic PDE problems?
Any comment would be appreciated.

[euan]

Quoting - minos67paran.com

Hi, I am trying to use the Pardiso solver for the nonstationary Navier-Stokes equations.
The only problem is that the system memory is increasing constantly.

I tried MKL_FreeBuffers(), but it was not working.
I am curious to know if the problem could be fixed.
Is anyone who use Pardiso solver to the dynamic PDE problems?
Any comment would be appreciated.

Have you called it with the phase parameter set to -1? Think you need to do that to free the memory

[minos67paran_com]

Quoting - euanw@smartmt.com

Have you called it with the phase parameter set to -1? Think you need to do that to free the memory

Yes, I did.
The following is a part of code calling pardiso.
I modified pardiso example code in MKL.
(I found that the c-code example usefortran array convention(starts index 1 for *ia, *ja).)
If there are somthings missing, please make some comment.


//----------------------------------------------------------

Sparse : custom structure for a sparse matrix
SparseToPardisoRowNum( mat, ia ); SparseToPardisoMat( mat, ja, a ); : custom function which I believe there
is no problem...

//-------------------------------------------------------------


int PARDISO_SOLVER_ICC(int nn, Sparse *mat, double *bb, double *xx )
{
/* Matrix data. */
MKL_INT n = ((MKL_INT) nn);
MKL_INT i, MKL_temp;
MKL_INT *ia, *ja;
double *a;

ia = (MKL_INT *)malloc( (n+1)*sizeof(MKL_INT) );
MKL_temp = SparseToPardisoRowNum( mat, ia );
// printf("check point -5. %dn", MKL_temp);

ja = (MKL_INT *)malloc( ia*sizeof(MKL_INT) );
a = (double *)malloc( ia*sizeof(double) );
MKL_temp = SparseToPardisoMat( mat, ja, a );
printf("Non-zeros of Pardiso mat = %dn", MKL_temp );
// printf("Non-zeros of Pardiso mat = %dn", SparseToPardisoMat( mat, ja, a ) );

// printf ("n");
// for (i = 0; i < n+1; i++) printf("ia[%d] = %dn", i, ia );
// printf ("n");
// for (i = 0; i < ia; i++) printf("ja[%d] = %dn", i, ja );
// printf ("n");
// for (i = 0; i < ia; i++) printf("a[%d] = %fn", i, a );

MKL_INT nnz = ia;
MKL_INT mtype = 11; /* Real unsymmetric matrix */

/* RHS and solution vectors. */
double *b, *x;
MKL_INT nrhs = 1; /* Number of right hand sides. */

b = (double *)malloc( n*sizeof(double) );
x = (double *)malloc( n*sizeof(double) );

/* 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];
int pt[64];

/* Pardiso control parameters. */
MKL_INT iparm[64];
MKL_INT maxfct, mnum, phase, error, msglvl;

/* Auxiliary variables. */

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 OMP_NUM_THREADS */
iparm[2] = 1;
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] = 1; /* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */
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 = 1; /* 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("check point 2. n" );

/* -------------------------------------------------------------------- */
/* .. 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 ...n ");
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
phase = 33;
// iparm[7] = 1; /* Max numbers of iterative refinement steps. */
iparm[7] = 2; /* Max numbers of iterative refinement steps. */

/* Set right hand side. */
for (i = 0; i < n; i++) {
b = bb;
}
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++) {
xx = x;
// printf("n x [%d] = %24.16e", i, x );
}
printf ("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);
return 0;
printf("n PARDISO_SOLVER done !");

MKL_FreeBuffers();

free(ia);
free(ja);

free(a);
free(b);
free(x);

}

[euan · Accepted solution]

hi,

i think you're doing the "return 0" at the end of your program before freeing your memory. I would have thought the compiler would have warned you of this, but then again I'm used to C# which tends to be more friendly

Quoting - minos67paran.com

Yes, I did.
The following is a part of code calling pardiso.
I modified pardiso example code in MKL.
(I found that the c-code example usefortran array convention(starts index 1 for *ia, *ja).)
If there are somthings missing, please make some comment.


//----------------------------------------------------------

Sparse : custom structure for a sparse matrix
SparseToPardisoRowNum( mat, ia ); SparseToPardisoMat( mat, ja, a ); : custom function which I believe there
is no problem...

//-------------------------------------------------------------


int PARDISO_SOLVER_ICC(int nn, Sparse *mat, double *bb, double *xx )
{
/* Matrix data. */
MKL_INT n = ((MKL_INT) nn);
MKL_INT i, MKL_temp;
MKL_INT *ia, *ja;
double *a;

ia = (MKL_INT *)malloc( (n+1)*sizeof(MKL_INT) );
MKL_temp = SparseToPardisoRowNum( mat, ia );
// printf("check point -5. %dn", MKL_temp);

ja = (MKL_INT *)malloc( ia*sizeof(MKL_INT) );
a = (double *)malloc( ia*sizeof(double) );
MKL_temp = SparseToPardisoMat( mat, ja, a );
printf("Non-zeros of Pardiso mat = %dn", MKL_temp );
// printf("Non-zeros of Pardiso mat = %dn", SparseToPardisoMat( mat, ja, a ) );

// printf ("n");
// for (i = 0; i < n+1; i++) printf("ia[%d] = %dn", i, ia );
// printf ("n");
// for (i = 0; i < ia; i++) printf("ja[%d] = %dn", i, ja );
// printf ("n");
// for (i = 0; i < ia; i++) printf("a[%d] = %fn", i, a );

MKL_INT nnz = ia;
MKL_INT mtype = 11; /* Real unsymmetric matrix */

/* RHS and solution vectors. */
double *b, *x;
MKL_INT nrhs = 1; /* Number of right hand sides. */

b = (double *)malloc( n*sizeof(double) );
x = (double *)malloc( n*sizeof(double) );

/* 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];
int pt[64];

/* Pardiso control parameters. */
MKL_INT iparm[64];
MKL_INT maxfct, mnum, phase, error, msglvl;

/* Auxiliary variables. */

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 OMP_NUM_THREADS */
iparm[2] = 1;
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] = 1; /* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */
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 = 1; /* 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("check point 2. n" );

/* -------------------------------------------------------------------- */
/* .. 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 ...n ");
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
phase = 33;
// iparm[7] = 1; /* Max numbers of iterative refinement steps. */
iparm[7] = 2; /* Max numbers of iterative refinement steps. */

/* Set right hand side. */
for (i = 0; i < n; i++) {
b = bb;
}
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++) {
xx = x;
// printf("n x [%d] = %24.16e", i, x );
}
printf ("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);
return 0;
printf("n PARDISO_SOLVER done !");

MKL_FreeBuffers();

free(ia);
free(ja);

free(a);
free(b);
free(x);

}

[minos67paran_com]

Quoting - euanw@smartmt.com

Thank you euanw!

After changing the function from int type to void type, the trouble is resolved.