Actual source code: mkl_pardiso.c

  1: #include <../src/mat/impls/aij/seq/aij.h>
  2: #include <../src/mat/impls/sbaij/seq/sbaij.h>
  3: #include <../src/mat/impls/dense/seq/dense.h>

  5: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
  6: #define MKL_ILP64
  7: #endif
  8: #include <mkl_pardiso.h>

 10: PETSC_EXTERN void PetscSetMKL_PARDISOThreads(int);

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

 28: #define IPARM_SIZE 64

 30: #if defined(PETSC_USE_64BIT_INDICES)
 31:  #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
 32:   #define INT_TYPE long long int
 33:   #define MKL_PARDISO pardiso
 34:   #define MKL_PARDISO_INIT pardisoinit
 35:  #else
 36:   /* this is the case where the MKL BLAS/LAPACK are 32 bit integers but the 64 bit integer version of
 37:      of Pardiso code is used; hence the need for the 64 below*/
 38:   #define INT_TYPE long long int
 39:   #define MKL_PARDISO pardiso_64
 40:   #define MKL_PARDISO_INIT pardiso_64init
 41: void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm [])
 42: {
 43:   int iparm_copy[IPARM_SIZE], mtype_copy, i;

 45:   mtype_copy = *mtype;
 46:   pardisoinit(pt, &mtype_copy, iparm_copy);
 47:   for (i=0; i<IPARM_SIZE; i++) iparm[i] = iparm_copy[i];
 48: }
 49:  #endif
 50: #else
 51:  #define INT_TYPE int
 52:  #define MKL_PARDISO pardiso
 53:  #define MKL_PARDISO_INIT pardisoinit
 54: #endif

 56: /*
 57:  *  Internal data structure.
 58:  *  For more information check mkl_pardiso manual.
 59:  */
 60: typedef struct {

 62:   /* Configuration vector*/
 63:   INT_TYPE     iparm[IPARM_SIZE];

 65:   /*
 66:    * Internal mkl_pardiso memory location.
 67:    * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
 68:    */
 69:   void         *pt[IPARM_SIZE];

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

 74:   /* Matrix structure*/
 75:   void         *a;
 76:   INT_TYPE     *ia, *ja;

 78:   /* Number of non-zero elements*/
 79:   INT_TYPE     nz;

 81:   /* Row permutaton vector*/
 82:   INT_TYPE     *perm;

 84:   /* Define if matrix preserves sparse structure.*/
 85:   MatStructure matstruc;

 87:   PetscBool    needsym;
 88:   PetscBool    freeaij;

 90:   /* Schur complement */
 91:   PetscScalar  *schur;
 92:   PetscInt     schur_size;
 93:   PetscInt     *schur_idxs;
 94:   PetscScalar  *schur_work;
 95:   PetscBLASInt schur_work_size;
 96:   PetscBool    solve_interior;

 98:   /* True if mkl_pardiso function have been used.*/
 99:   PetscBool CleanUp;

101:   /* Conversion to a format suitable for MKL */
102:   PetscErrorCode (*Convert)(Mat, PetscBool, MatReuse, PetscBool*, INT_TYPE*, INT_TYPE**, INT_TYPE**, PetscScalar**);
103: } Mat_MKL_PARDISO;

105: PetscErrorCode MatMKLPardiso_Convert_seqsbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
106: {
107:   Mat_SeqSBAIJ   *aa = (Mat_SeqSBAIJ*)A->data;
108:   PetscInt       bs  = A->rmap->bs,i;

112:   if (!sym) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
113:   *v      = aa->a;
114:   if (bs == 1) { /* already in the correct format */
115:     /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
116:     *r    = (INT_TYPE*)aa->i;
117:     *c    = (INT_TYPE*)aa->j;
118:     *nnz  = (INT_TYPE)aa->nz;
119:     *free = PETSC_FALSE;
120:   } else if (reuse == MAT_INITIAL_MATRIX) {
121:     PetscInt m = A->rmap->n,nz = aa->nz;
122:     PetscInt *row,*col;
123:     PetscMalloc2(m+1,&row,nz,&col);
124:     for (i=0; i<m+1; i++) {
125:       row[i] = aa->i[i]+1;
126:     }
127:     for (i=0; i<nz; i++) {
128:       col[i] = aa->j[i]+1;
129:     }
130:     *r    = (INT_TYPE*)row;
131:     *c    = (INT_TYPE*)col;
132:     *nnz  = (INT_TYPE)nz;
133:     *free = PETSC_TRUE;
134:   }
135:   return(0);
136: }

138: PetscErrorCode MatMKLPardiso_Convert_seqbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
139: {
140:   Mat_SeqBAIJ    *aa = (Mat_SeqBAIJ*)A->data;
141:   PetscInt       bs  = A->rmap->bs,i;

145:   if (!sym) {
146:     *v      = aa->a;
147:     if (bs == 1) { /* already in the correct format */
148:       /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
149:       *r    = (INT_TYPE*)aa->i;
150:       *c    = (INT_TYPE*)aa->j;
151:       *nnz  = (INT_TYPE)aa->nz;
152:       *free = PETSC_FALSE;
153:       return(0);
154:     } else if (reuse == MAT_INITIAL_MATRIX) {
155:       PetscInt m = A->rmap->n,nz = aa->nz;
156:       PetscInt *row,*col;
157:       PetscMalloc2(m+1,&row,nz,&col);
158:       for (i=0; i<m+1; i++) {
159:         row[i] = aa->i[i]+1;
160:       }
161:       for (i=0; i<nz; i++) {
162:         col[i] = aa->j[i]+1;
163:       }
164:       *r    = (INT_TYPE*)row;
165:       *c    = (INT_TYPE*)col;
166:       *nnz  = (INT_TYPE)nz;
167:     }
168:     *free = PETSC_TRUE;
169:   } else {
170:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
171:   }
172:   return(0);
173: }

175: PetscErrorCode MatMKLPardiso_Convert_seqaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
176: {
177:   Mat_SeqAIJ     *aa = (Mat_SeqAIJ*)A->data;
178:   PetscScalar    *aav;

182:   MatSeqAIJGetArrayRead(A,(const PetscScalar**)&aav);
183:   if (!sym) { /* already in the correct format */
184:     *v    = aav;
185:     *r    = (INT_TYPE*)aa->i;
186:     *c    = (INT_TYPE*)aa->j;
187:     *nnz  = (INT_TYPE)aa->nz;
188:     *free = PETSC_FALSE;
189:   } else if (reuse == MAT_INITIAL_MATRIX) { /* need to get the triangular part */
190:     PetscScalar *vals,*vv;
191:     PetscInt    *row,*col,*jj;
192:     PetscInt    m = A->rmap->n,nz,i;

194:     nz = 0;
195:     for (i=0; i<m; i++) nz += aa->i[i+1] - aa->diag[i];
196:     PetscMalloc2(m+1,&row,nz,&col);
197:     PetscMalloc1(nz,&vals);
198:     jj = col;
199:     vv = vals;

201:     row[0] = 0;
202:     for (i=0; i<m; i++) {
203:       PetscInt    *aj = aa->j + aa->diag[i];
204:       PetscScalar *av = aav + aa->diag[i];
205:       PetscInt    rl  = aa->i[i+1] - aa->diag[i],j;

207:       for (j=0; j<rl; j++) {
208:         *jj = *aj; jj++; aj++;
209:         *vv = *av; vv++; av++;
210:       }
211:       row[i+1] = row[i] + rl;
212:     }
213:     *v    = vals;
214:     *r    = (INT_TYPE*)row;
215:     *c    = (INT_TYPE*)col;
216:     *nnz  = (INT_TYPE)nz;
217:     *free = PETSC_TRUE;
218:   } else {
219:     PetscScalar *vv;
220:     PetscInt    m = A->rmap->n,i;

222:     vv = *v;
223:     for (i=0; i<m; i++) {
224:       PetscScalar *av = aav + aa->diag[i];
225:       PetscInt    rl  = aa->i[i+1] - aa->diag[i],j;
226:       for (j=0; j<rl; j++) {
227:         *vv = *av; vv++; av++;
228:       }
229:     }
230:     *free = PETSC_TRUE;
231:   }
232:   MatSeqAIJRestoreArrayRead(A,(const PetscScalar**)&aav);
233:   return(0);
234: }

236: static PetscErrorCode MatMKLPardisoSolveSchur_Private(Mat F, PetscScalar *B, PetscScalar *X)
237: {
238:   Mat_MKL_PARDISO      *mpardiso = (Mat_MKL_PARDISO*)F->data;
239:   Mat                  S,Xmat,Bmat;
240:   MatFactorSchurStatus schurstatus;
241:   PetscErrorCode       ierr;

244:   MatFactorGetSchurComplement(F,&S,&schurstatus);
245:   if (X == B && schurstatus == MAT_FACTOR_SCHUR_INVERTED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"X and B cannot point to the same address");
246:   MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,B,&Bmat);
247:   MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,X,&Xmat);
248:   MatSetType(Bmat,((PetscObject)S)->type_name);
249:   MatSetType(Xmat,((PetscObject)S)->type_name);
250: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
251:   MatBindToCPU(Xmat,S->boundtocpu);
252:   MatBindToCPU(Bmat,S->boundtocpu);
253: #endif

255: #if defined(PETSC_USE_COMPLEX)
256:   if (mpardiso->iparm[12-1] == 1) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Hermitian solve not implemented yet");
257: #endif

259:   switch (schurstatus) {
260:   case MAT_FACTOR_SCHUR_FACTORED:
261:     if (!mpardiso->iparm[12-1]) {
262:       MatMatSolve(S,Bmat,Xmat);
263:     } else { /* transpose solve */
264:       MatMatSolveTranspose(S,Bmat,Xmat);
265:     }
266:     break;
267:   case MAT_FACTOR_SCHUR_INVERTED:
268:     MatProductCreateWithMat(S,Bmat,NULL,Xmat);
269:     if (!mpardiso->iparm[12-1]) {
270:       MatProductSetType(Xmat,MATPRODUCT_AB);
271:     } else { /* transpose solve */
272:       MatProductSetType(Xmat,MATPRODUCT_AtB);
273:     }
274:     MatProductSetFromOptions(Xmat);
275:     MatProductSymbolic(Xmat);
276:     MatProductNumeric(Xmat);
277:     MatProductClear(Xmat);
278:     break;
279:   default:
280:     SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status);
281:     break;
282:   }
283:   MatFactorRestoreSchurComplement(F,&S,schurstatus);
284:   MatDestroy(&Bmat);
285:   MatDestroy(&Xmat);
286:   return(0);
287: }

289: PetscErrorCode MatFactorSetSchurIS_MKL_PARDISO(Mat F, IS is)
290: {
291:   Mat_MKL_PARDISO   *mpardiso = (Mat_MKL_PARDISO*)F->data;
292:   const PetscScalar *arr;
293:   const PetscInt    *idxs;
294:   PetscInt          size,i;
295:   PetscMPIInt       csize;
296:   PetscBool         sorted;
297:   PetscErrorCode    ierr;

300:   MPI_Comm_size(PetscObjectComm((PetscObject)F),&csize);
301:   if (csize > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MKL_PARDISO parallel Schur complements not yet supported from PETSc");
302:   ISSorted(is,&sorted);
303:   if (!sorted) {
304:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS for MKL_PARDISO Schur complements needs to be sorted");
305:   }
306:   ISGetLocalSize(is,&size);
307:   PetscFree(mpardiso->schur_work);
308:   PetscBLASIntCast(PetscMax(mpardiso->n,2*size),&mpardiso->schur_work_size);
309:   PetscMalloc1(mpardiso->schur_work_size,&mpardiso->schur_work);
310:   MatDestroy(&F->schur);
311:   MatCreateSeqDense(PETSC_COMM_SELF,size,size,NULL,&F->schur);
312:   MatDenseGetArrayRead(F->schur,&arr);
313:   mpardiso->schur      = (PetscScalar*)arr;
314:   mpardiso->schur_size = size;
315:   MatDenseRestoreArrayRead(F->schur,&arr);
316:   if (mpardiso->mtype == 2) {
317:     MatSetOption(F->schur,MAT_SPD,PETSC_TRUE);
318:   }

320:   PetscFree(mpardiso->schur_idxs);
321:   PetscMalloc1(size,&mpardiso->schur_idxs);
322:   PetscArrayzero(mpardiso->perm,mpardiso->n);
323:   ISGetIndices(is,&idxs);
324:   PetscArraycpy(mpardiso->schur_idxs,idxs,size);
325:   for (i=0;i<size;i++) mpardiso->perm[idxs[i]] = 1;
326:   ISRestoreIndices(is,&idxs);
327:   if (size) { /* turn on Schur switch if the set of indices is not empty */
328:     mpardiso->iparm[36-1] = 2;
329:   }
330:   return(0);
331: }

333: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
334: {
335:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
336:   PetscErrorCode  ierr;

339:   if (mat_mkl_pardiso->CleanUp) {
340:     mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

342:     MKL_PARDISO (mat_mkl_pardiso->pt,
343:       &mat_mkl_pardiso->maxfct,
344:       &mat_mkl_pardiso->mnum,
345:       &mat_mkl_pardiso->mtype,
346:       &mat_mkl_pardiso->phase,
347:       &mat_mkl_pardiso->n,
348:       NULL,
349:       NULL,
350:       NULL,
351:       NULL,
352:       &mat_mkl_pardiso->nrhs,
353:       mat_mkl_pardiso->iparm,
354:       &mat_mkl_pardiso->msglvl,
355:       NULL,
356:       NULL,
357:       &mat_mkl_pardiso->err);
358:   }
359:   PetscFree(mat_mkl_pardiso->perm);
360:   PetscFree(mat_mkl_pardiso->schur_work);
361:   PetscFree(mat_mkl_pardiso->schur_idxs);
362:   if (mat_mkl_pardiso->freeaij) {
363:     PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
364:     if (mat_mkl_pardiso->iparm[34] == 1) {
365:       PetscFree(mat_mkl_pardiso->a);
366:     }
367:   }
368:   PetscFree(A->data);

370:   /* clear composed functions */
371:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
372:   PetscObjectComposeFunction((PetscObject)A,"MatFactorSetSchurIS_C",NULL);
373:   PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);
374:   return(0);
375: }

377: static PetscErrorCode MatMKLPardisoScatterSchur_Private(Mat_MKL_PARDISO *mpardiso, PetscScalar *whole, PetscScalar *schur, PetscBool reduce)
378: {
380:   if (reduce) { /* data given for the whole matrix */
381:     PetscInt i,m=0,p=0;
382:     for (i=0;i<mpardiso->nrhs;i++) {
383:       PetscInt j;
384:       for (j=0;j<mpardiso->schur_size;j++) {
385:         schur[p+j] = whole[m+mpardiso->schur_idxs[j]];
386:       }
387:       m += mpardiso->n;
388:       p += mpardiso->schur_size;
389:     }
390:   } else { /* from Schur to whole */
391:     PetscInt i,m=0,p=0;
392:     for (i=0;i<mpardiso->nrhs;i++) {
393:       PetscInt j;
394:       for (j=0;j<mpardiso->schur_size;j++) {
395:         whole[m+mpardiso->schur_idxs[j]] = schur[p+j];
396:       }
397:       m += mpardiso->n;
398:       p += mpardiso->schur_size;
399:     }
400:   }
401:   return(0);
402: }

404: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x)
405: {
406:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
407:   PetscErrorCode    ierr;
408:   PetscScalar       *xarray;
409:   const PetscScalar *barray;

412:   mat_mkl_pardiso->nrhs = 1;
413:   VecGetArrayWrite(x,&xarray);
414:   VecGetArrayRead(b,&barray);

416:   if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
417:   else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;

419:   if (barray == xarray) { /* if the two vectors share the same memory */
420:     PetscScalar *work;
421:     if (!mat_mkl_pardiso->schur_work) {
422:       PetscMalloc1(mat_mkl_pardiso->n,&work);
423:     } else {
424:       work = mat_mkl_pardiso->schur_work;
425:     }
426:     mat_mkl_pardiso->iparm[6-1] = 1;
427:     MKL_PARDISO (mat_mkl_pardiso->pt,
428:       &mat_mkl_pardiso->maxfct,
429:       &mat_mkl_pardiso->mnum,
430:       &mat_mkl_pardiso->mtype,
431:       &mat_mkl_pardiso->phase,
432:       &mat_mkl_pardiso->n,
433:       mat_mkl_pardiso->a,
434:       mat_mkl_pardiso->ia,
435:       mat_mkl_pardiso->ja,
436:       NULL,
437:       &mat_mkl_pardiso->nrhs,
438:       mat_mkl_pardiso->iparm,
439:       &mat_mkl_pardiso->msglvl,
440:       (void*)xarray,
441:       (void*)work,
442:       &mat_mkl_pardiso->err);
443:     if (!mat_mkl_pardiso->schur_work) {
444:       PetscFree(work);
445:     }
446:   } else {
447:     mat_mkl_pardiso->iparm[6-1] = 0;
448:     MKL_PARDISO (mat_mkl_pardiso->pt,
449:       &mat_mkl_pardiso->maxfct,
450:       &mat_mkl_pardiso->mnum,
451:       &mat_mkl_pardiso->mtype,
452:       &mat_mkl_pardiso->phase,
453:       &mat_mkl_pardiso->n,
454:       mat_mkl_pardiso->a,
455:       mat_mkl_pardiso->ia,
456:       mat_mkl_pardiso->ja,
457:       mat_mkl_pardiso->perm,
458:       &mat_mkl_pardiso->nrhs,
459:       mat_mkl_pardiso->iparm,
460:       &mat_mkl_pardiso->msglvl,
461:       (void*)barray,
462:       (void*)xarray,
463:       &mat_mkl_pardiso->err);
464:   }
465:   VecRestoreArrayRead(b,&barray);

467:   if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);

469:   if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
470:     if (!mat_mkl_pardiso->solve_interior) {
471:       PetscInt shift = mat_mkl_pardiso->schur_size;

473:       MatFactorFactorizeSchurComplement(A);
474:       /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
475:       if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;

477:       /* solve Schur complement */
478:       MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
479:       MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
480:       MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
481:     } else { /* if we are solving for the interior problem, any value in barray[schur] forward-substituted to xarray[schur] will be neglected */
482:       PetscInt i;
483:       for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
484:         xarray[mat_mkl_pardiso->schur_idxs[i]] = 0.;
485:       }
486:     }

488:     /* expansion phase */
489:     mat_mkl_pardiso->iparm[6-1] = 1;
490:     mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
491:     MKL_PARDISO (mat_mkl_pardiso->pt,
492:       &mat_mkl_pardiso->maxfct,
493:       &mat_mkl_pardiso->mnum,
494:       &mat_mkl_pardiso->mtype,
495:       &mat_mkl_pardiso->phase,
496:       &mat_mkl_pardiso->n,
497:       mat_mkl_pardiso->a,
498:       mat_mkl_pardiso->ia,
499:       mat_mkl_pardiso->ja,
500:       mat_mkl_pardiso->perm,
501:       &mat_mkl_pardiso->nrhs,
502:       mat_mkl_pardiso->iparm,
503:       &mat_mkl_pardiso->msglvl,
504:       (void*)xarray,
505:       (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
506:       &mat_mkl_pardiso->err);

508:     if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
509:     mat_mkl_pardiso->iparm[6-1] = 0;
510:   }
511:   VecRestoreArrayWrite(x,&xarray);
512:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
513:   return(0);
514: }

516: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x)
517: {
518:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
519:   PetscInt        oiparm12;
520:   PetscErrorCode  ierr;

523:   oiparm12 = mat_mkl_pardiso->iparm[12 - 1];
524:   mat_mkl_pardiso->iparm[12 - 1] = 2;
525:   MatSolve_MKL_PARDISO(A,b,x);
526:   mat_mkl_pardiso->iparm[12 - 1] = oiparm12;
527:   return(0);
528: }

530: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X)
531: {
532:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->data;
533:   PetscErrorCode    ierr;
534:   const PetscScalar *barray;
535:   PetscScalar       *xarray;
536:   PetscBool         flg;

539:   PetscObjectBaseTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
540:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
541:   if (X != B) {
542:     PetscObjectBaseTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
543:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");
544:   }

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

548:   if (mat_mkl_pardiso->nrhs > 0) {
549:     MatDenseGetArrayRead(B,&barray);
550:     MatDenseGetArrayWrite(X,&xarray);

552:     if (barray == xarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"B and X cannot share the same memory location");
553:     if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
554:     else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;

556:     MKL_PARDISO (mat_mkl_pardiso->pt,
557:       &mat_mkl_pardiso->maxfct,
558:       &mat_mkl_pardiso->mnum,
559:       &mat_mkl_pardiso->mtype,
560:       &mat_mkl_pardiso->phase,
561:       &mat_mkl_pardiso->n,
562:       mat_mkl_pardiso->a,
563:       mat_mkl_pardiso->ia,
564:       mat_mkl_pardiso->ja,
565:       mat_mkl_pardiso->perm,
566:       &mat_mkl_pardiso->nrhs,
567:       mat_mkl_pardiso->iparm,
568:       &mat_mkl_pardiso->msglvl,
569:       (void*)barray,
570:       (void*)xarray,
571:       &mat_mkl_pardiso->err);
572:     if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);

574:     MatDenseRestoreArrayRead(B,&barray);
575:     if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
576:       PetscScalar *o_schur_work = NULL;

578:       /* solve Schur complement */
579:       if (!mat_mkl_pardiso->solve_interior) {
580:         PetscInt shift = mat_mkl_pardiso->schur_size*mat_mkl_pardiso->nrhs,scale;
581:         PetscInt mem = mat_mkl_pardiso->n*mat_mkl_pardiso->nrhs;

583:         MatFactorFactorizeSchurComplement(A);
584:         /* allocate extra memory if it is needed */
585:         scale = 1;
586:         if (A->schur_status == MAT_FACTOR_SCHUR_INVERTED) scale = 2;
587:         mem *= scale;
588:         if (mem > mat_mkl_pardiso->schur_work_size) {
589:           o_schur_work = mat_mkl_pardiso->schur_work;
590:           PetscMalloc1(mem,&mat_mkl_pardiso->schur_work);
591:         }
592:         /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
593:         if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;
594:         MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
595:         MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
596:         MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
597:       } else { /* if we are solving for the interior problem, any value in barray[schur,n] forward-substituted to xarray[schur,n] will be neglected */
598:         PetscInt i,n,m=0;
599:         for (n=0;n<mat_mkl_pardiso->nrhs;n++) {
600:           for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
601:             xarray[mat_mkl_pardiso->schur_idxs[i]+m] = 0.;
602:           }
603:           m += mat_mkl_pardiso->n;
604:         }
605:       }

607:       /* expansion phase */
608:       mat_mkl_pardiso->iparm[6-1] = 1;
609:       mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
610:       MKL_PARDISO (mat_mkl_pardiso->pt,
611:         &mat_mkl_pardiso->maxfct,
612:         &mat_mkl_pardiso->mnum,
613:         &mat_mkl_pardiso->mtype,
614:         &mat_mkl_pardiso->phase,
615:         &mat_mkl_pardiso->n,
616:         mat_mkl_pardiso->a,
617:         mat_mkl_pardiso->ia,
618:         mat_mkl_pardiso->ja,
619:         mat_mkl_pardiso->perm,
620:         &mat_mkl_pardiso->nrhs,
621:         mat_mkl_pardiso->iparm,
622:         &mat_mkl_pardiso->msglvl,
623:         (void*)xarray,
624:         (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
625:         &mat_mkl_pardiso->err);
626:       if (o_schur_work) { /* restore original schur_work (minimal size) */
627:         PetscFree(mat_mkl_pardiso->schur_work);
628:         mat_mkl_pardiso->schur_work = o_schur_work;
629:       }
630:       if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
631:       mat_mkl_pardiso->iparm[6-1] = 0;
632:     }
633:     MatDenseRestoreArrayWrite(X,&xarray);
634:   }
635:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
636:   return(0);
637: }

639: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info)
640: {
641:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->data;
642:   PetscErrorCode  ierr;

645:   mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
646:   (*mat_mkl_pardiso->Convert)(A,mat_mkl_pardiso->needsym,MAT_REUSE_MATRIX,&mat_mkl_pardiso->freeaij,&mat_mkl_pardiso->nz,&mat_mkl_pardiso->ia,&mat_mkl_pardiso->ja,(PetscScalar**)&mat_mkl_pardiso->a);

648:   mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
649:   MKL_PARDISO (mat_mkl_pardiso->pt,
650:     &mat_mkl_pardiso->maxfct,
651:     &mat_mkl_pardiso->mnum,
652:     &mat_mkl_pardiso->mtype,
653:     &mat_mkl_pardiso->phase,
654:     &mat_mkl_pardiso->n,
655:     mat_mkl_pardiso->a,
656:     mat_mkl_pardiso->ia,
657:     mat_mkl_pardiso->ja,
658:     mat_mkl_pardiso->perm,
659:     &mat_mkl_pardiso->nrhs,
660:     mat_mkl_pardiso->iparm,
661:     &mat_mkl_pardiso->msglvl,
662:     NULL,
663:     (void*)mat_mkl_pardiso->schur,
664:     &mat_mkl_pardiso->err);
665:   if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);

667:   /* report flops */
668:   if (mat_mkl_pardiso->iparm[18] > 0) {
669:     PetscLogFlops(PetscPowRealInt(10.,6)*mat_mkl_pardiso->iparm[18]);
670:   }

672:   if (F->schur) { /* schur output from pardiso is in row major format */
673: #if defined(PETSC_HAVE_CUDA)
674:     F->schur->offloadmask = PETSC_OFFLOAD_CPU;
675: #endif
676:     MatFactorRestoreSchurComplement(F,NULL,MAT_FACTOR_SCHUR_UNFACTORED);
677:     MatTranspose(F->schur,MAT_INPLACE_MATRIX,&F->schur);
678:   }
679:   mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
680:   mat_mkl_pardiso->CleanUp  = PETSC_TRUE;
681:   return(0);
682: }

684: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A)
685: {
686:   Mat_MKL_PARDISO     *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
687:   PetscErrorCode      ierr;
688:   PetscInt            icntl,bs,threads=1;
689:   PetscBool           flg;

692:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");

694:   PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use within PARDISO","None",threads,&threads,&flg);
695:   if (flg) PetscSetMKL_PARDISOThreads((int)threads);

697:   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);
698:   if (flg) mat_mkl_pardiso->maxfct = icntl;

700:   PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);
701:   if (flg) mat_mkl_pardiso->mnum = icntl;

703:   PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);
704:   if (flg) mat_mkl_pardiso->msglvl = icntl;

706:   PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);
707:   if (flg) {
708:     void *pt[IPARM_SIZE];
709:     mat_mkl_pardiso->mtype = icntl;
710:     icntl = mat_mkl_pardiso->iparm[34];
711:     bs = mat_mkl_pardiso->iparm[36];
712:     MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
713: #if defined(PETSC_USE_REAL_SINGLE)
714:     mat_mkl_pardiso->iparm[27] = 1;
715: #else
716:     mat_mkl_pardiso->iparm[27] = 0;
717: #endif
718:     mat_mkl_pardiso->iparm[34] = icntl;
719:     mat_mkl_pardiso->iparm[36] = bs;
720:   }

722:   PetscOptionsInt("-mat_mkl_pardiso_1","Use default values (if 0)","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);
723:   if (flg) mat_mkl_pardiso->iparm[0] = icntl;

725:   PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);
726:   if (flg) mat_mkl_pardiso->iparm[1] = icntl;

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

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

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

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

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

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

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

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

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

755:   PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations (0 to disable)","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);
756:   if (flg) mat_mkl_pardiso->iparm[18] = icntl;

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

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

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

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

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

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

776:   PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
777:   if (flg) mat_mkl_pardiso->iparm[59] = icntl;
778:   PetscOptionsEnd();
779:   return(0);
780: }

782: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
783: {
785:   PetscInt       i,bs;
786:   PetscBool      match;

789:   for (i=0; i<IPARM_SIZE; i++) mat_mkl_pardiso->iparm[i] = 0;
790:   for (i=0; i<IPARM_SIZE; i++) mat_mkl_pardiso->pt[i] = 0;
791: #if defined(PETSC_USE_REAL_SINGLE)
792:   mat_mkl_pardiso->iparm[27] = 1;
793: #else
794:   mat_mkl_pardiso->iparm[27] = 0;
795: #endif
796:   /* Default options for both sym and unsym */
797:   mat_mkl_pardiso->iparm[ 0] =  1; /* Solver default parameters overriden with provided by iparm */
798:   mat_mkl_pardiso->iparm[ 1] =  2; /* Metis reordering */
799:   mat_mkl_pardiso->iparm[ 5] =  0; /* Write solution into x */
800:   mat_mkl_pardiso->iparm[ 7] =  0; /* Max number of iterative refinement steps */
801:   mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
802:   mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
803: #if 0
804:   mat_mkl_pardiso->iparm[23] =  1; /* Parallel factorization control*/
805: #endif
806:   PetscObjectTypeCompareAny((PetscObject)A,&match,MATSEQBAIJ,MATSEQSBAIJ,"");
807:   MatGetBlockSize(A,&bs);
808:   if (!match || bs == 1) {
809:     mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
810:     mat_mkl_pardiso->n         = A->rmap->N;
811:   } else {
812:     mat_mkl_pardiso->iparm[34] = 0; /* Cluster Sparse Solver use Fortran-style indexing for ia and ja arrays */
813:     mat_mkl_pardiso->iparm[36] = bs;
814:     mat_mkl_pardiso->n         = A->rmap->N/bs;
815:   }
816:   mat_mkl_pardiso->iparm[39] =  0; /* Input: matrix/rhs/solution stored on rank-0 */

818:   mat_mkl_pardiso->CleanUp   = PETSC_FALSE;
819:   mat_mkl_pardiso->maxfct    = 1; /* Maximum number of numerical factorizations. */
820:   mat_mkl_pardiso->mnum      = 1; /* Which factorization to use. */
821:   mat_mkl_pardiso->msglvl    = 0; /* 0: do not print 1: Print statistical information in file */
822:   mat_mkl_pardiso->phase     = -1;
823:   mat_mkl_pardiso->err       = 0;

825:   mat_mkl_pardiso->nrhs      = 1;
826:   mat_mkl_pardiso->err       = 0;
827:   mat_mkl_pardiso->phase     = -1;

829:   if (ftype == MAT_FACTOR_LU) {
830:     mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
831:     mat_mkl_pardiso->iparm[10] =  1; /* Use nonsymmetric permutation and scaling MPS */
832:     mat_mkl_pardiso->iparm[12] =  1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
833:   } else {
834:     mat_mkl_pardiso->iparm[ 9] = 8; /* Perturb the pivot elements with 1E-8 */
835:     mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
836:     mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
837: #if defined(PETSC_USE_DEBUG)
838:     mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
839: #endif
840:   }
841:   PetscCalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
842:   mat_mkl_pardiso->schur_size = 0;
843:   return(0);
844: }

846: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
847: {
848:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
849:   PetscErrorCode  ierr;

852:   mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
853:   PetscSetMKL_PARDISOFromOptions(F,A);
854:   /* throw away any previously computed structure */
855:   if (mat_mkl_pardiso->freeaij) {
856:     PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
857:     if (mat_mkl_pardiso->iparm[34] == 1) {
858:       PetscFree(mat_mkl_pardiso->a);
859:     }
860:   }
861:   (*mat_mkl_pardiso->Convert)(A,mat_mkl_pardiso->needsym,MAT_INITIAL_MATRIX,&mat_mkl_pardiso->freeaij,&mat_mkl_pardiso->nz,&mat_mkl_pardiso->ia,&mat_mkl_pardiso->ja,(PetscScalar**)&mat_mkl_pardiso->a);
862:   if (mat_mkl_pardiso->iparm[34] == 1) mat_mkl_pardiso->n = A->rmap->N;
863:   else mat_mkl_pardiso->n = A->rmap->N/A->rmap->bs;

865:   mat_mkl_pardiso->phase = JOB_ANALYSIS;

867:   /* reset flops counting if requested */
868:   if (mat_mkl_pardiso->iparm[18]) mat_mkl_pardiso->iparm[18] = -1;

870:   MKL_PARDISO (mat_mkl_pardiso->pt,
871:     &mat_mkl_pardiso->maxfct,
872:     &mat_mkl_pardiso->mnum,
873:     &mat_mkl_pardiso->mtype,
874:     &mat_mkl_pardiso->phase,
875:     &mat_mkl_pardiso->n,
876:     mat_mkl_pardiso->a,
877:     mat_mkl_pardiso->ia,
878:     mat_mkl_pardiso->ja,
879:     mat_mkl_pardiso->perm,
880:     &mat_mkl_pardiso->nrhs,
881:     mat_mkl_pardiso->iparm,
882:     &mat_mkl_pardiso->msglvl,
883:     NULL,
884:     NULL,
885:     &mat_mkl_pardiso->err);
886:   if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);

888:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;

890:   if (F->factortype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
891:   else F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO;

893:   F->ops->solve           = MatSolve_MKL_PARDISO;
894:   F->ops->solvetranspose  = MatSolveTranspose_MKL_PARDISO;
895:   F->ops->matsolve        = MatMatSolve_MKL_PARDISO;
896:   return(0);
897: }

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

904:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
905:   return(0);
906: }

908: #if !defined(PETSC_USE_COMPLEX)
909: PetscErrorCode MatGetInertia_MKL_PARDISO(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
910: {
911:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)F->data;

914:   if (nneg) *nneg = mat_mkl_pardiso->iparm[22];
915:   if (npos) *npos = mat_mkl_pardiso->iparm[21];
916:   if (nzero) *nzero = F->rmap->N - (mat_mkl_pardiso->iparm[22] + mat_mkl_pardiso->iparm[21]);
917:   return(0);
918: }
919: #endif

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

926:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
927: #if defined(PETSC_USE_COMPLEX)
928:   F->ops->getinertia = NULL;
929: #else
930:   F->ops->getinertia = MatGetInertia_MKL_PARDISO;
931: #endif
932:   return(0);
933: }

935: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
936: {
937:   PetscErrorCode    ierr;
938:   PetscBool         iascii;
939:   PetscViewerFormat format;
940:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
941:   PetscInt          i;

944:   if (A->ops->solve != MatSolve_MKL_PARDISO) return(0);

946:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
947:   if (iascii) {
948:     PetscViewerGetFormat(viewer,&format);
949:     if (format == PETSC_VIEWER_ASCII_INFO) {
950:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
951:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase:             %d \n",mat_mkl_pardiso->phase);
952:       for (i=1; i<=64; i++) {
953:         PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]:     %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
954:       }
955:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct:     %d \n", mat_mkl_pardiso->maxfct);
956:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum:     %d \n", mat_mkl_pardiso->mnum);
957:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype:     %d \n", mat_mkl_pardiso->mtype);
958:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n:     %d \n", mat_mkl_pardiso->n);
959:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs:     %d \n", mat_mkl_pardiso->nrhs);
960:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl:     %d \n", mat_mkl_pardiso->msglvl);
961:     }
962:   }
963:   return(0);
964: }

966: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
967: {
968:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)A->data;

971:   info->block_size        = 1.0;
972:   info->nz_used           = mat_mkl_pardiso->iparm[17];
973:   info->nz_allocated      = mat_mkl_pardiso->iparm[17];
974:   info->nz_unneeded       = 0.0;
975:   info->assemblies        = 0.0;
976:   info->mallocs           = 0.0;
977:   info->memory            = 0.0;
978:   info->fill_ratio_given  = 0;
979:   info->fill_ratio_needed = 0;
980:   info->factor_mallocs    = 0;
981:   return(0);
982: }

984: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
985: {
986:   PetscInt        backup,bs;
987:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;

990:   if (icntl <= 64) {
991:     mat_mkl_pardiso->iparm[icntl - 1] = ival;
992:   } else {
993:     if (icntl == 65) PetscSetMKL_PARDISOThreads(ival);
994:     else if (icntl == 66) mat_mkl_pardiso->maxfct = ival;
995:     else if (icntl == 67) mat_mkl_pardiso->mnum = ival;
996:     else if (icntl == 68) mat_mkl_pardiso->msglvl = ival;
997:     else if (icntl == 69) {
998:       void *pt[IPARM_SIZE];
999:       backup = mat_mkl_pardiso->iparm[34];
1000:       bs = mat_mkl_pardiso->iparm[36];
1001:       mat_mkl_pardiso->mtype = ival;
1002:       MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
1003: #if defined(PETSC_USE_REAL_SINGLE)
1004:       mat_mkl_pardiso->iparm[27] = 1;
1005: #else
1006:       mat_mkl_pardiso->iparm[27] = 0;
1007: #endif
1008:       mat_mkl_pardiso->iparm[34] = backup;
1009:       mat_mkl_pardiso->iparm[36] = bs;
1010:     } else if (icntl==70) mat_mkl_pardiso->solve_interior = (PetscBool)!!ival;
1011:   }
1012:   return(0);
1013: }

1015: /*@
1016:   MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters

1018:    Logically Collective on Mat

1020:    Input Parameters:
1021: +  F - the factored matrix obtained by calling MatGetFactor()
1022: .  icntl - index of Mkl_Pardiso parameter
1023: -  ival - value of Mkl_Pardiso parameter

1025:   Options Database:
1026: .   -mat_mkl_pardiso_<icntl> <ival>

1028:    Level: beginner

1030:    References:
1031: .      Mkl_Pardiso Users' Guide

1033: .seealso: MatGetFactor()
1034: @*/
1035: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
1036: {

1040:   PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
1041:   return(0);
1042: }

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

1048:   Works with MATSEQAIJ matrices

1050:   Use -pc_type lu -pc_factor_mat_solver_type mkl_pardiso to use this direct solver

1052:   Options Database Keys:
1053: + -mat_mkl_pardiso_65 - Number of threads to use within MKL_PARDISO
1054: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
1055: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
1056: . -mat_mkl_pardiso_68 - Message level information
1057: . -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
1058: . -mat_mkl_pardiso_1  - Use default values
1059: . -mat_mkl_pardiso_2  - Fill-in reducing ordering for the input matrix
1060: . -mat_mkl_pardiso_4  - Preconditioned CGS/CG
1061: . -mat_mkl_pardiso_5  - User permutation
1062: . -mat_mkl_pardiso_6  - Write solution on x
1063: . -mat_mkl_pardiso_8  - Iterative refinement step
1064: . -mat_mkl_pardiso_10 - Pivoting perturbation
1065: . -mat_mkl_pardiso_11 - Scaling vectors
1066: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
1067: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
1068: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
1069: . -mat_mkl_pardiso_19 - Report number of floating point operations
1070: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
1071: . -mat_mkl_pardiso_24 - Parallel factorization control
1072: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
1073: . -mat_mkl_pardiso_27 - Matrix checker
1074: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
1075: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
1076: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode

1078:   Level: beginner

1080:   For more information please check  mkl_pardiso manual

1082: .seealso: PCFactorSetMatSolverType(), MatSolverType

1084: M*/
1085: static PetscErrorCode MatFactorGetSolverType_mkl_pardiso(Mat A, MatSolverType *type)
1086: {
1088:   *type = MATSOLVERMKL_PARDISO;
1089:   return(0);
1090: }

1092: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
1093: {
1094:   Mat             B;
1095:   PetscErrorCode  ierr;
1096:   Mat_MKL_PARDISO *mat_mkl_pardiso;
1097:   PetscBool       isSeqAIJ,isSeqBAIJ,isSeqSBAIJ;

1100:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1101:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1102:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1103:   MatCreate(PetscObjectComm((PetscObject)A),&B);
1104:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1105:   PetscStrallocpy("mkl_pardiso",&((PetscObject)B)->type_name);
1106:   MatSetUp(B);

1108:   PetscNewLog(B,&mat_mkl_pardiso);
1109:   B->data = mat_mkl_pardiso;

1111:   MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
1112:   if (ftype == MAT_FACTOR_LU) {
1113:     B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
1114:     B->factortype            = MAT_FACTOR_LU;
1115:     mat_mkl_pardiso->needsym = PETSC_FALSE;
1116:     if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1117:     else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1118:     else if (isSeqSBAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU factor with SEQSBAIJ format! Use MAT_FACTOR_CHOLESKY instead");
1119:     else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU with %s format",((PetscObject)A)->type_name);
1120: #if defined(PETSC_USE_COMPLEX)
1121:     mat_mkl_pardiso->mtype = 13;
1122: #else
1123:     mat_mkl_pardiso->mtype = 11;
1124: #endif
1125:   } else {
1126:     B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
1127:     B->factortype                  = MAT_FACTOR_CHOLESKY;
1128:     if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1129:     else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1130:     else if (isSeqSBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqsbaij;
1131:     else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with %s format",((PetscObject)A)->type_name);

1133:     mat_mkl_pardiso->needsym = PETSC_TRUE;
1134: #if !defined(PETSC_USE_COMPLEX)
1135:     if (A->spd_set && A->spd) mat_mkl_pardiso->mtype = 2;
1136:     else                      mat_mkl_pardiso->mtype = -2;
1137: #else
1138:     mat_mkl_pardiso->mtype = 6;
1139:     if (A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with Hermitian matrices! Use MAT_FACTOR_LU instead");
1140: #endif
1141:   }
1142:   B->ops->destroy = MatDestroy_MKL_PARDISO;
1143:   B->ops->view    = MatView_MKL_PARDISO;
1144:   B->ops->getinfo = MatGetInfo_MKL_PARDISO;
1145:   B->factortype   = ftype;
1146:   B->assembled    = PETSC_TRUE;

1148:   PetscFree(B->solvertype);
1149:   PetscStrallocpy(MATSOLVERMKL_PARDISO,&B->solvertype);

1151:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_pardiso);
1152:   PetscObjectComposeFunction((PetscObject)B,"MatFactorSetSchurIS_C",MatFactorSetSchurIS_MKL_PARDISO);
1153:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);

1155:   *F = B;
1156:   return(0);
1157: }

1159: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void)
1160: {

1164:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1165:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1166:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQBAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1167:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1168:   return(0);
1169: }