Here are some of the codes
void CPoissonSolver::solve_lineareqn_mkl_symm()
{
int jid, kid, nid, nn0, icount;
nn0 = matVals_.size();
int matrix_size = ndof_;
int *ia = new int[ matrix_size+1]; //next row start index
int *ja = new int[nn0];
double *a = new double[nn0];
double *b = new double[matrix_size];
double *x = new double[matrix_size]; /* RHS and solution vectors. */
icount = 0;
ia[0] = 1;
for(map<int,set<int>>::iterator it = matIds_.begin(); it != matIds_.end(); ++it )
{
jid = it->first;
for(set<int>::iterator sit = it->second.begin(); sit != it->second.end(); ++sit)
{
kid = *sit;
nid = jid * ndof_ + kid;
a[icount] = matVals_[nid];
ja[icount] = kid+1;
icount++;
}
ia[jid+1] = icount+1; //jav.size()+1;
// cout<<jid<<" :: "<<icount<<" "<<ia[jid+1]<<endl;
}
cout<<"symmetric matrix!!!:";
cout<<" non zeros : "<<nn0<<" "<<icount<<endl;
for (int i=0;i < matrix_size;++i)
b[i] = rhs_[i];
/* cout<<matrix_size<<" :: ";
for (int i=0;i <= matrix_size;++i)
{
cout<<ia[i]<<" ";
}
cout<<endl;
for (int i=0; i < nn0;++i)
{
cout<<a[i]<<" "<<ja[i]<<endl;
}*/
int nrhs = 1; /* Number of right hand sides. */
int mtype = -2; /* Real symmetric matrix */
void *pt[64]; /* Internal solver memory pointer pt, */
int iparm[64]; /* Pardiso control parameters. */
int maxfct, mnum, phase, error, msglvl;
int i;
double ddum; /* Double dummy *//* Auxiliary variables. */
int idum; /* Integer dummy. */
for (i = 0; i < 64; i++)
iparm[i] = 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;//omp_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] = 16; /* Default logical fortran unit number for output */
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; /* Not in use */
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; /* Don't 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[i] = 0;
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
phase = 11; //only reordering and symbolic factorization
PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
&matrix_size, 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]);