Actual source code: mkl_pardiso.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64)
  2: #define MKL_ILP64
  3: #endif

  5: #include <../src/mat/impls/aij/seq/aij.h>    /*I "petscmat.h" I*/
  6: #include <../src/mat/impls/dense/seq/dense.h>

  8: #include <stdio.h>
  9: #include <stdlib.h>
 10: #include <math.h>
 11: #include <mkl.h>

 13: /*
 14:  *  Possible mkl_pardiso phases that controls the execution of the solver.
 15:  *  For more information check mkl_pardiso manual.
 16:  */
 17: #define JOB_ANALYSIS 11
 18: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
 19: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
 20: #define JOB_NUMERICAL_FACTORIZATION 22
 21: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
 22: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
 23: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
 24: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
 25: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
 26: #define JOB_RELEASE_OF_LU_MEMORY 0
 27: #define JOB_RELEASE_OF_ALL_MEMORY -1

 29: #define IPARM_SIZE 64

 31: #if defined(PETSC_USE_64BIT_INDICES)
 32:  #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64)
 33:   /* sizeof(MKL_INT) == sizeof(long long int) if ilp64*/
 34:   #define INT_TYPE long long int
 35:   #define MKL_PARDISO pardiso
 36:   #define MKL_PARDISO_INIT pardisoinit
 37:  #else
 38:   #define INT_TYPE long long int
 39:   #define MKL_PARDISO pardiso_64
 40:   #define MKL_PARDISO_INIT pardiso_64init
 41:  #endif
 42: #else
 43:  #define INT_TYPE int
 44:  #define MKL_PARDISO pardiso
 45:  #define MKL_PARDISO_INIT pardisoinit
 46: #endif


 49: /*
 50:  *  Internal data structure.
 51:  *  For more information check mkl_pardiso manual.
 52:  */
 53: typedef struct {

 55:   /* Configuration vector*/
 56:   INT_TYPE     iparm[IPARM_SIZE];

 58:   /*
 59:    * Internal mkl_pardiso memory location.
 60:    * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
 61:    */
 62:   void         *pt[IPARM_SIZE];

 64:   /* Basic mkl_pardiso info*/
 65:   INT_TYPE     phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;

 67:   /* Matrix structure*/
 68:   void         *a;
 69:   INT_TYPE     *ia, *ja;

 71:   /* Number of non-zero elements*/
 72:   INT_TYPE     nz;

 74:   /* Row permutaton vector*/
 75:   INT_TYPE     *perm;

 77:   /* Define if matrix preserves sparse structure.*/
 78:   MatStructure matstruc;

 80:   /* True if mkl_pardiso function have been used.*/
 81:   PetscBool CleanUp;
 82: } Mat_MKL_PARDISO;


 85: void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm [])
 86: {
 87:   int iparm_copy[IPARM_SIZE], mtype_copy, i;
 88: 
 89:   mtype_copy = *mtype;
 90:   pardisoinit(pt, &mtype_copy, iparm_copy);
 91:   for(i = 0; i < IPARM_SIZE; i++){
 92:     iparm[i] = iparm_copy[i];
 93:   }
 94: }


 97: /*
 98:  * Copy the elements of matrix A.
 99:  * Input:
100:  *   - Mat A: MATSEQAIJ matrix
101:  *   - int shift: matrix index.
102:  *     - 0 for c representation
103:  *     - 1 for fortran representation
104:  *   - MatReuse reuse:
105:  *     - MAT_INITIAL_MATRIX: Create a new aij representation
106:  *     - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
107:  * Output:
108:  *   - int *nnz: Number of nonzero-elements.
109:  *   - int **r pointer to i index
110:  *   - int **c pointer to j elements
111:  *   - MATRIXTYPE **v: Non-zero elements
112:  */
115: PetscErrorCode MatCopy_MKL_PARDISO(Mat A, MatReuse reuse, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, void **v)
116: {
117:   Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data;

120:   *v=aa->a;
121:   if (reuse == MAT_INITIAL_MATRIX) {
122:     *r   = (INT_TYPE*)aa->i;
123:     *c   = (INT_TYPE*)aa->j;
124:     *nnz = aa->nz;
125:   }
126:   return(0);
127: }

129: /*
130:  * Free memory for Mat_MKL_PARDISO structure and pointers to objects.
131:  */
134: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
135: {
136:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
137:   PetscErrorCode  ierr;

140:   /* Terminate instance, deallocate memories */
141:   if (mat_mkl_pardiso->CleanUp) {
142:     mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

144:     MKL_PARDISO (mat_mkl_pardiso->pt,
145:       &mat_mkl_pardiso->maxfct,
146:       &mat_mkl_pardiso->mnum,
147:       &mat_mkl_pardiso->mtype,
148:       &mat_mkl_pardiso->phase,
149:       &mat_mkl_pardiso->n,
150:       NULL,
151:       NULL,
152:       NULL,
153:       mat_mkl_pardiso->perm,
154:       &mat_mkl_pardiso->nrhs,
155:       mat_mkl_pardiso->iparm,
156:       &mat_mkl_pardiso->msglvl,
157:       NULL,
158:       NULL,
159:       &mat_mkl_pardiso->err);
160:   }
161:   PetscFree(mat_mkl_pardiso->perm);
162:   PetscFree(A->spptr);

164:   /* clear composed functions */
165:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
166:   PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);

168:   MatDestroy_SeqAIJ(A);
169:   return(0);
170: }

172: /*
173:  * Computes Ax = b
174:  */
177: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x)
178: {
179:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
180:   PetscErrorCode    ierr;
181:   PetscScalar       *xarray;
182:   const PetscScalar *barray;

185:   mat_mkl_pardiso->nrhs = 1;
186:   VecGetArray(x,&xarray);
187:   VecGetArrayRead(b,&barray);

189:   /* solve phase */
190:   /*-------------*/
191:   mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
192:   MKL_PARDISO (mat_mkl_pardiso->pt,
193:     &mat_mkl_pardiso->maxfct,
194:     &mat_mkl_pardiso->mnum,
195:     &mat_mkl_pardiso->mtype,
196:     &mat_mkl_pardiso->phase,
197:     &mat_mkl_pardiso->n,
198:     mat_mkl_pardiso->a,
199:     mat_mkl_pardiso->ia,
200:     mat_mkl_pardiso->ja,
201:     mat_mkl_pardiso->perm,
202:     &mat_mkl_pardiso->nrhs,
203:     mat_mkl_pardiso->iparm,
204:     &mat_mkl_pardiso->msglvl,
205:     (void*)barray,
206:     (void*)xarray,
207:     &mat_mkl_pardiso->err);

209:   if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err);
210:   VecRestoreArray(x,&xarray);
211:   VecRestoreArrayRead(b,&barray);
212:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
213:   return(0);
214: }


219: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x)
220: {
221:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
222:   PetscErrorCode  ierr;

225: #if defined(PETSC_USE_COMPLEX)
226:   mat_mkl_pardiso->iparm[12 - 1] = 1;
227: #else
228:   mat_mkl_pardiso->iparm[12 - 1] = 2;
229: #endif
230:   MatSolve_MKL_PARDISO(A,b,x);
231:   mat_mkl_pardiso->iparm[12 - 1] = 0;
232:   return(0);
233: }


238: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X)
239: {
240:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
241:   PetscErrorCode    ierr;
242:   PetscScalar       *barray, *xarray;
243:   PetscBool         flg;

246:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
247:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
248:   PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
249:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");

251:   MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);

253:   if(mat_mkl_pardiso->nrhs > 0){
254:     MatDenseGetArray(B,&barray);
255:     MatDenseGetArray(X,&xarray);

257:     /* solve phase */
258:     /*-------------*/
259:     mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
260:     MKL_PARDISO (mat_mkl_pardiso->pt,
261:       &mat_mkl_pardiso->maxfct,
262:       &mat_mkl_pardiso->mnum,
263:       &mat_mkl_pardiso->mtype,
264:       &mat_mkl_pardiso->phase,
265:       &mat_mkl_pardiso->n,
266:       mat_mkl_pardiso->a,
267:       mat_mkl_pardiso->ia,
268:       mat_mkl_pardiso->ja,
269:       mat_mkl_pardiso->perm,
270:       &mat_mkl_pardiso->nrhs,
271:       mat_mkl_pardiso->iparm,
272:       &mat_mkl_pardiso->msglvl,
273:       (void*)barray,
274:       (void*)xarray,
275:       &mat_mkl_pardiso->err);
276:     if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err);
277:   }
278:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
279:   return(0);
280: }

282: /*
283:  * LU Decomposition
284:  */
287: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info)
288: {
289:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->spptr;
290:   PetscErrorCode  ierr;

292:   /* numerical factorization phase */
293:   /*-------------------------------*/
295:   mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
296:   MatCopy_MKL_PARDISO(A, MAT_REUSE_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);

298:   /* numerical factorization phase */
299:   /*-------------------------------*/
300:   mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
301:   MKL_PARDISO (mat_mkl_pardiso->pt,
302:     &mat_mkl_pardiso->maxfct,
303:     &mat_mkl_pardiso->mnum,
304:     &mat_mkl_pardiso->mtype,
305:     &mat_mkl_pardiso->phase,
306:     &mat_mkl_pardiso->n,
307:     mat_mkl_pardiso->a,
308:     mat_mkl_pardiso->ia,
309:     mat_mkl_pardiso->ja,
310:     mat_mkl_pardiso->perm,
311:     &mat_mkl_pardiso->nrhs,
312:     mat_mkl_pardiso->iparm,
313:     &mat_mkl_pardiso->msglvl,
314:     NULL,
315:     NULL,
316:     &mat_mkl_pardiso->err);
317:   if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err);

319:   mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
320:   mat_mkl_pardiso->CleanUp  = PETSC_TRUE;
321:   return(0);
322: }

324: /* Sets mkl_pardiso options from the options database */
327: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A)
328: {
329:   Mat_MKL_PARDISO     *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
330:   PetscErrorCode      ierr;
331:   PetscInt            icntl;
332:   PetscBool           flg;
333:   int                 pt[IPARM_SIZE], threads = 1;

336:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");
337:   PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use","None",threads,&threads,&flg);
338:   if (flg) mkl_set_num_threads(threads);

340:   PetscOptionsInt("-mat_mkl_pardiso_66","Maximum number of factors with identical sparsity structure that must be kept in memory at the same time","None",mat_mkl_pardiso->maxfct,&icntl,&flg);
341:   if (flg) mat_mkl_pardiso->maxfct = icntl;

343:   PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);
344:   if (flg) mat_mkl_pardiso->mnum = icntl;
345: 
346:   PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);
347:   if (flg) mat_mkl_pardiso->msglvl = icntl;

349:   PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);
350:   if(flg){
351:    mat_mkl_pardiso->mtype = icntl;
352:    MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
353: #if defined(PETSC_USE_REAL_SINGLE)
354:     mat_mkl_pardiso->iparm[27] = 1;
355: #else
356:     mat_mkl_pardiso->iparm[27] = 0;
357: #endif
358:     mat_mkl_pardiso->iparm[34] = 1;
359:   }
360:   PetscOptionsInt("-mat_mkl_pardiso_1","Use default values","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);

362:   if(flg && icntl != 0){
363:     PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);
364:     if (flg) mat_mkl_pardiso->iparm[1] = icntl;

366:     PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);
367:     if (flg) mat_mkl_pardiso->iparm[3] = icntl;

369:     PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);
370:     if (flg) mat_mkl_pardiso->iparm[4] = icntl;

372:     PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);
373:     if (flg) mat_mkl_pardiso->iparm[5] = icntl;

375:     PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);
376:     if (flg) mat_mkl_pardiso->iparm[7] = icntl;

378:     PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);
379:     if (flg) mat_mkl_pardiso->iparm[9] = icntl;

381:     PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);
382:     if (flg) mat_mkl_pardiso->iparm[10] = icntl;

384:     PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);
385:     if (flg) mat_mkl_pardiso->iparm[11] = icntl;

387:     PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);
388:     if (flg) mat_mkl_pardiso->iparm[12] = icntl;

390:     PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);
391:     if (flg) mat_mkl_pardiso->iparm[17] = icntl;

393:     PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);
394:     if (flg) mat_mkl_pardiso->iparm[18] = icntl;

396:     PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);
397:     if (flg) mat_mkl_pardiso->iparm[20] = icntl;

399:     PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);
400:     if (flg) mat_mkl_pardiso->iparm[23] = icntl;

402:     PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);
403:     if (flg) mat_mkl_pardiso->iparm[24] = icntl;

405:     PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);
406:     if (flg) mat_mkl_pardiso->iparm[26] = icntl;

408:     PetscOptionsInt("-mat_mkl_pardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_pardiso->iparm[30],&icntl,&flg);
409:     if (flg) mat_mkl_pardiso->iparm[30] = icntl;

411:     PetscOptionsInt("-mat_mkl_pardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_pardiso->iparm[33],&icntl,&flg);
412:     if (flg) mat_mkl_pardiso->iparm[33] = icntl;

414:     PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
415:     if (flg) mat_mkl_pardiso->iparm[59] = icntl;
416:   }
417:   PetscOptionsEnd();
418:   return(0);
419: }

423: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
424: {
426:   PetscInt       i;

429:   for ( i = 0; i < IPARM_SIZE; i++ ){
430:     mat_mkl_pardiso->iparm[i] = 0;
431:   }

433:   for ( i = 0; i < IPARM_SIZE; i++ ){
434:     mat_mkl_pardiso->pt[i] = 0;
435:   }
436: 
437:   /* Default options for both sym and unsym */
438:   mat_mkl_pardiso->iparm[ 0] =  1; /* Solver default parameters overriden with provided by iparm */
439:   mat_mkl_pardiso->iparm[ 1] =  2; /* Metis reordering */
440:   mat_mkl_pardiso->iparm[ 5] =  0; /* Write solution into x */
441:   mat_mkl_pardiso->iparm[ 7] =  2; /* Max number of iterative refinement steps */
442:   mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
443:   mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
444: #if 0
445:   mat_mkl_pardiso->iparm[23] =  1; /* Parallel factorization control*/
446: #endif
447:   mat_mkl_pardiso->iparm[34] =  1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
448:   mat_mkl_pardiso->iparm[39] =  0; /* Input: matrix/rhs/solution stored on master */
449: 
450:   mat_mkl_pardiso->CleanUp   = PETSC_FALSE;
451:   mat_mkl_pardiso->maxfct    = 1; /* Maximum number of numerical factorizations. */
452:   mat_mkl_pardiso->mnum      = 1; /* Which factorization to use. */
453:   mat_mkl_pardiso->msglvl    = 0; /* 0: do not print 1: Print statistical information in file */
454:   mat_mkl_pardiso->phase     = -1;
455:   mat_mkl_pardiso->err       = 0;
456: 
457:   mat_mkl_pardiso->n         = A->rmap->N;
458:   mat_mkl_pardiso->nrhs      = 1;
459:   mat_mkl_pardiso->err       = 0;
460:   mat_mkl_pardiso->phase     = -1;
461: 
462:   if(ftype == MAT_FACTOR_LU){
463:     /* Default type for non-sym */
464: #if defined(PETSC_USE_COMPLEX)
465:     mat_mkl_pardiso->mtype     = 13;
466: #else
467:     mat_mkl_pardiso->mtype     = 11;
468: #endif

470:     mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
471:     mat_mkl_pardiso->iparm[10] =  1; /* Use nonsymmetric permutation and scaling MPS */
472:     mat_mkl_pardiso->iparm[12] =  1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */

474:   } else {
475:     /* Default type for sym */
476: #if defined(PETSC_USE_COMPLEX)
477:     mat_mkl_pardiso ->mtype    = 3;
478: #else
479:     mat_mkl_pardiso ->mtype    = -2;
480: #endif
481:     mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
482:     mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
483:     mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
484: /*    mat_mkl_pardiso->iparm[20] =  1; */ /* Apply 1x1 and 2x2 Bunch-Kaufman pivoting during the factorization process */
485: #if defined(PETSC_USE_DEBUG)
486:     mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
487: #endif
488:   }
489:   PetscMalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
490:   for(i = 0; i < A->rmap->N; i++){
491:     mat_mkl_pardiso->perm[i] = 0;
492:   }
493:   return(0);
494: }

496: /*
497:  * Symbolic decomposition. Mkl_Pardiso analysis phase.
498:  */
501: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
502: {
503:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
504:   PetscErrorCode  ierr;

507:   mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;

509:   /* Set MKL_PARDISO options from the options database */
510:   PetscSetMKL_PARDISOFromOptions(F,A);

512:   MatCopy_MKL_PARDISO(A, MAT_INITIAL_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);
513:   mat_mkl_pardiso->n = A->rmap->N;

515:   /* analysis phase */
516:   /*----------------*/
517:   mat_mkl_pardiso->phase = JOB_ANALYSIS;

519:   MKL_PARDISO (mat_mkl_pardiso->pt,
520:     &mat_mkl_pardiso->maxfct,
521:     &mat_mkl_pardiso->mnum,
522:     &mat_mkl_pardiso->mtype,
523:     &mat_mkl_pardiso->phase,
524:     &mat_mkl_pardiso->n,
525:     mat_mkl_pardiso->a,
526:     mat_mkl_pardiso->ia,
527:     mat_mkl_pardiso->ja,
528:     mat_mkl_pardiso->perm,
529:     &mat_mkl_pardiso->nrhs,
530:     mat_mkl_pardiso->iparm,
531:     &mat_mkl_pardiso->msglvl,
532:     NULL,
533:     NULL,
534:     &mat_mkl_pardiso->err);
535:   if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d\n. Please check manual",mat_mkl_pardiso->err);

537:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;

539:   if(F->factortype == MAT_FACTOR_LU){
540:     F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
541:   } else {
542:     F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO;
543:   }
544:   F->ops->solve           = MatSolve_MKL_PARDISO;
545:   F->ops->solvetranspose  = MatSolveTranspose_MKL_PARDISO;
546:   F->ops->matsolve        = MatMatSolve_MKL_PARDISO;
547:   return(0);
548: }

552: PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
553: {

557:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
558:   return(0);
559: }

563: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,const MatFactorInfo *info)
564: {

568:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
569:   return(0);
570: }

574: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
575: {
576:   PetscErrorCode    ierr;
577:   PetscBool         iascii;
578:   PetscViewerFormat format;
579:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
580:   PetscInt          i;

583:   /* check if matrix is mkl_pardiso type */
584:   if (A->ops->solve != MatSolve_MKL_PARDISO) return(0);

586:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
587:   if (iascii) {
588:     PetscViewerGetFormat(viewer,&format);
589:     if (format == PETSC_VIEWER_ASCII_INFO) {
590:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
591:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase:             %d \n",mat_mkl_pardiso->phase);
592:       for(i = 1; i <= 64; i++){
593:         PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]:     %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
594:       }
595:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct:     %d \n", mat_mkl_pardiso->maxfct);
596:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum:     %d \n", mat_mkl_pardiso->mnum);
597:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype:     %d \n", mat_mkl_pardiso->mtype);
598:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n:     %d \n", mat_mkl_pardiso->n);
599:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs:     %d \n", mat_mkl_pardiso->nrhs);
600:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl:     %d \n", mat_mkl_pardiso->msglvl);
601:     }
602:   }
603:   return(0);
604: }


609: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
610: {
611:   Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->spptr;

614:   info->block_size        = 1.0;
615:   info->nz_allocated      = mat_mkl_pardiso->nz + 0.0;
616:   info->nz_unneeded       = 0.0;
617:   info->assemblies        = 0.0;
618:   info->mallocs           = 0.0;
619:   info->memory            = 0.0;
620:   info->fill_ratio_given  = 0;
621:   info->fill_ratio_needed = 0;
622:   info->factor_mallocs    = 0;
623:   return(0);
624: }

628: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
629: {
630:   Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->spptr;

633:   if(icntl <= 64){
634:     mat_mkl_pardiso->iparm[icntl - 1] = ival;
635:   } else {
636:     if(icntl == 65)
637:       mkl_set_num_threads((int)ival);
638:     else if(icntl == 66)
639:       mat_mkl_pardiso->maxfct = ival;
640:     else if(icntl == 67)
641:       mat_mkl_pardiso->mnum = ival;
642:     else if(icntl == 68)
643:       mat_mkl_pardiso->msglvl = ival;
644:     else if(icntl == 69){
645:       int pt[IPARM_SIZE];
646:       mat_mkl_pardiso->mtype = ival;
647:       MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
648: #if defined(PETSC_USE_REAL_SINGLE)
649:       mat_mkl_pardiso->iparm[27] = 1;
650: #else
651:       mat_mkl_pardiso->iparm[27] = 0;
652: #endif
653:       mat_mkl_pardiso->iparm[34] = 1;
654:     }
655:   }
656:   return(0);
657: }

661: /*@
662:   MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters

664:    Logically Collective on Mat

666:    Input Parameters:
667: +  F - the factored matrix obtained by calling MatGetFactor()
668: .  icntl - index of Mkl_Pardiso parameter
669: -  ival - value of Mkl_Pardiso parameter

671:   Options Database:
672: .   -mat_mkl_pardiso_<icntl> <ival>

674:    Level: beginner

676:    References: Mkl_Pardiso Users' Guide

678: .seealso: MatGetFactor()
679: @*/
680: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
681: {

685:   PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
686:   return(0);
687: }

689: /*MC
690:   MATSOLVERMKL_PARDISO -  A matrix type providing direct solvers (LU) for
691:   sequential matrices via the external package MKL_PARDISO.

693:   Works with MATSEQAIJ matrices

695:   Use -pc_type lu -pc_factor_mat_solver_package mkl_pardiso to us this direct solver

697:   Options Database Keys:
698: + -mat_mkl_pardiso_65 - Number of thrads to use
699: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
700: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
701: . -mat_mkl_pardiso_68 - Message level information
702: . -mat_mkl_pardiso_69 - Defines the matrix type. IMPORTANT: When you set this flag, iparm parameters are going to be set to the default ones for the matrix type
703: . -mat_mkl_pardiso_1 - Use default values
704: . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix
705: . -mat_mkl_pardiso_4 - Preconditioned CGS/CG
706: . -mat_mkl_pardiso_5 - User permutation
707: . -mat_mkl_pardiso_6 - Write solution on x
708: . -mat_mkl_pardiso_8 - Iterative refinement step
709: . -mat_mkl_pardiso_10 - Pivoting perturbation
710: . -mat_mkl_pardiso_11 - Scaling vectors
711: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
712: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
713: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
714: . -mat_mkl_pardiso_19 - Report number of floating point operations
715: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
716: . -mat_mkl_pardiso_24 - Parallel factorization control
717: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
718: . -mat_mkl_pardiso_27 - Matrix checker
719: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
720: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
721: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode

723:   Level: beginner

725:   For more information please check  mkl_pardiso manual

727: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage

729: M*/
732: static PetscErrorCode MatFactorGetSolverPackage_mkl_pardiso(Mat A, const MatSolverPackage *type)
733: {
735:   *type = MATSOLVERMKL_PARDISO;
736:   return(0);
737: }

739: /* MatGetFactor for Seq sbAIJ matrices */
742: PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
743: {
744:   Mat             B;
745:   PetscErrorCode  ierr;
746:   Mat_MKL_PARDISO *mat_mkl_pardiso;
747:   PetscBool       isSeqSBAIJ;
748:   PetscInt        bs;

751:   /* Create the factorization matrix */
752:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
753:   MatCreate(PetscObjectComm((PetscObject)A),&B);
754:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
755:   MatSetType(B,((PetscObject)A)->type_name);
756:   MatSeqSBAIJSetPreallocation(B,1,0,NULL);
757:   MatGetBlockSize(A,&bs);

759:   if(bs != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQSBAIJ with block size other than 1 is not supported by Pardiso");
760:   if(ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_CHOLESKY.");
761: 
762:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
763:   B->factortype                  = MAT_FACTOR_CHOLESKY;
764:   B->ops->destroy                = MatDestroy_MKL_PARDISO;
765:   B->ops->view                   = MatView_MKL_PARDISO;
766:   B->factortype                  = ftype;
767:   B->ops->getinfo                = MatGetInfo_MKL_PARDISO;
768:   B->assembled                   = PETSC_TRUE;           /* required by -ksp_view */

770:   PetscNewLog(B,&mat_mkl_pardiso);
771:   B->spptr = mat_mkl_pardiso;
772:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);
773:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
774:   MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
775:   *F = B;
776:   return(0);
777: }

779: /* MatGetFactor for Seq AIJ matrices */
782: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
783: {
784:   Mat             B;
785:   PetscErrorCode  ierr;
786:   Mat_MKL_PARDISO *mat_mkl_pardiso;
787:   PetscBool       isSeqAIJ;

790:   /* Create the factorization matrix */
791:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
792:   MatCreate(PetscObjectComm((PetscObject)A),&B);
793:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
794:   MatSetType(B,((PetscObject)A)->type_name);
795:   MatSeqAIJSetPreallocation(B,0,NULL);

797:   if(ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_LU.");

799:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
800:   B->factortype            = MAT_FACTOR_LU;
801:   B->ops->destroy          = MatDestroy_MKL_PARDISO;
802:   B->ops->view             = MatView_MKL_PARDISO;
803:   B->factortype            = ftype;
804:   B->ops->getinfo          = MatGetInfo_MKL_PARDISO;
805:   B->assembled             = PETSC_TRUE;           /* required by -ksp_view */

807:   PetscNewLog(B,&mat_mkl_pardiso);
808:   B->spptr = mat_mkl_pardiso;
809:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);
810:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
811:   MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);

813:   *F = B;
814:   return(0);
815: }

819: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MKL_Pardiso(void)
820: {

824:   MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,   MAT_FACTOR_LU,      MatGetFactor_aij_mkl_pardiso  );
825:   MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mkl_pardiso);
826:   return(0);
827: }