Actual source code: mkl_pardiso.c

petsc-master 2020-07-09
Report Typos and Errors
  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


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

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

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

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

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

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

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

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

 88:   PetscBool    needsym;
 89:   PetscBool    freeaij;

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

336: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
337: {
338:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
339:   PetscErrorCode  ierr;

342:   if (mat_mkl_pardiso->CleanUp) {
343:     mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

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

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

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

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

415:   mat_mkl_pardiso->nrhs = 1;
416:   VecGetArray(x,&xarray);
417:   VecGetArrayRead(b,&barray);

419:   if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
420:   else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;

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

470:   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);

472:   if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
473:     PetscInt shift = mat_mkl_pardiso->schur_size;

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

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

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

510:     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);
511:     mat_mkl_pardiso->iparm[6-1] = 0;
512:   }
513:   VecRestoreArray(x,&xarray);
514:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
515:   return(0);
516: }

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

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

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

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

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

550:   if (mat_mkl_pardiso->nrhs > 0) {
551:     MatDenseGetArrayRead(B,&barray);
552:     MatDenseGetArray(X,&xarray);

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

558:     MKL_PARDISO (mat_mkl_pardiso->pt,
559:       &mat_mkl_pardiso->maxfct,
560:       &mat_mkl_pardiso->mnum,
561:       &mat_mkl_pardiso->mtype,
562:       &mat_mkl_pardiso->phase,
563:       &mat_mkl_pardiso->n,
564:       mat_mkl_pardiso->a,
565:       mat_mkl_pardiso->ia,
566:       mat_mkl_pardiso->ja,
567:       mat_mkl_pardiso->perm,
568:       &mat_mkl_pardiso->nrhs,
569:       mat_mkl_pardiso->iparm,
570:       &mat_mkl_pardiso->msglvl,
571:       (void*)barray,
572:       (void*)xarray,
573:       &mat_mkl_pardiso->err);
574:     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);

576:     MatDenseRestoreArrayRead(B,&barray);
577:     if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
578:       PetscScalar *o_schur_work = NULL;
579:       PetscInt    shift = mat_mkl_pardiso->schur_size*mat_mkl_pardiso->nrhs,scale;
580:       PetscInt    mem = mat_mkl_pardiso->n*mat_mkl_pardiso->nrhs;

582:       /* allocate extra memory if it is needed */
583:       scale = 1;
584:       if (A->schur_status == MAT_FACTOR_SCHUR_INVERTED) scale = 2;

586:       mem *= scale;
587:       if (mem > mat_mkl_pardiso->schur_work_size) {
588:         o_schur_work = mat_mkl_pardiso->schur_work;
589:         PetscMalloc1(mem,&mat_mkl_pardiso->schur_work);
590:       }

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

595:       /* solve Schur complement */
596:       if (!mat_mkl_pardiso->solve_interior) {
597:         MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
598:         MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
599:         MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
600:       } else { /* if we are solving for the interior problem, any value in barray[schur,n] forward-substituted to xarray[schur,n] will be neglected */
601:         PetscInt i,n,m=0;
602:         for (n=0;n<mat_mkl_pardiso->nrhs;n++) {
603:           for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
604:             xarray[mat_mkl_pardiso->schur_idxs[i]+m] = 0.;
605:           }
606:           m += mat_mkl_pardiso->n;
607:         }
608:       }

610:       /* expansion phase */
611:       mat_mkl_pardiso->iparm[6-1] = 1;
612:       mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
613:       MKL_PARDISO (mat_mkl_pardiso->pt,
614:         &mat_mkl_pardiso->maxfct,
615:         &mat_mkl_pardiso->mnum,
616:         &mat_mkl_pardiso->mtype,
617:         &mat_mkl_pardiso->phase,
618:         &mat_mkl_pardiso->n,
619:         mat_mkl_pardiso->a,
620:         mat_mkl_pardiso->ia,
621:         mat_mkl_pardiso->ja,
622:         mat_mkl_pardiso->perm,
623:         &mat_mkl_pardiso->nrhs,
624:         mat_mkl_pardiso->iparm,
625:         &mat_mkl_pardiso->msglvl,
626:         (void*)xarray,
627:         (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
628:         &mat_mkl_pardiso->err);
629:       if (o_schur_work) { /* restore original schur_work (minimal size) */
630:         PetscFree(mat_mkl_pardiso->schur_work);
631:         mat_mkl_pardiso->schur_work = o_schur_work;
632:       }
633:       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);
634:       mat_mkl_pardiso->iparm[6-1] = 0;
635:     }
636:   }
637:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
638:   return(0);
639: }

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

647:   mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
648:   (*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);

650:   mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
651:   MKL_PARDISO (mat_mkl_pardiso->pt,
652:     &mat_mkl_pardiso->maxfct,
653:     &mat_mkl_pardiso->mnum,
654:     &mat_mkl_pardiso->mtype,
655:     &mat_mkl_pardiso->phase,
656:     &mat_mkl_pardiso->n,
657:     mat_mkl_pardiso->a,
658:     mat_mkl_pardiso->ia,
659:     mat_mkl_pardiso->ja,
660:     mat_mkl_pardiso->perm,
661:     &mat_mkl_pardiso->nrhs,
662:     mat_mkl_pardiso->iparm,
663:     &mat_mkl_pardiso->msglvl,
664:     NULL,
665:     (void*)mat_mkl_pardiso->schur,
666:     &mat_mkl_pardiso->err);
667:   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);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

784: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
785: {
787:   PetscInt       i,bs;
788:   PetscBool      match;

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

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

827:   mat_mkl_pardiso->nrhs      = 1;
828:   mat_mkl_pardiso->err       = 0;
829:   mat_mkl_pardiso->phase     = -1;

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

851: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
852: {
853:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
854:   PetscErrorCode  ierr;

857:   mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
858:   PetscSetMKL_PARDISOFromOptions(F,A);

860:   /* throw away any previously computed structure */
861:   if (mat_mkl_pardiso->freeaij) {
862:     PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
863:     if (mat_mkl_pardiso->iparm[34] == 1) {
864:       PetscFree(mat_mkl_pardiso->a);
865:     }
866:   }
867:   (*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);
868:   if (mat_mkl_pardiso->iparm[34] == 1) mat_mkl_pardiso->n = A->rmap->N;
869:   else mat_mkl_pardiso->n = A->rmap->N/A->rmap->bs;

871:   mat_mkl_pardiso->phase = JOB_ANALYSIS;

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

876:   MKL_PARDISO (mat_mkl_pardiso->pt,
877:     &mat_mkl_pardiso->maxfct,
878:     &mat_mkl_pardiso->mnum,
879:     &mat_mkl_pardiso->mtype,
880:     &mat_mkl_pardiso->phase,
881:     &mat_mkl_pardiso->n,
882:     mat_mkl_pardiso->a,
883:     mat_mkl_pardiso->ia,
884:     mat_mkl_pardiso->ja,
885:     mat_mkl_pardiso->perm,
886:     &mat_mkl_pardiso->nrhs,
887:     mat_mkl_pardiso->iparm,
888:     &mat_mkl_pardiso->msglvl,
889:     NULL,
890:     NULL,
891:     &mat_mkl_pardiso->err);
892:   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);

894:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;

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

899:   F->ops->solve           = MatSolve_MKL_PARDISO;
900:   F->ops->solvetranspose  = MatSolveTranspose_MKL_PARDISO;
901:   F->ops->matsolve        = MatMatSolve_MKL_PARDISO;
902:   return(0);
903: }

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

910:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
911:   return(0);
912: }

914: #if !defined(PETSC_USE_COMPLEX)
915: PetscErrorCode MatGetInertia_MKL_PARDISO(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
916: {
917:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)F->data;

920:   if (nneg) *nneg = mat_mkl_pardiso->iparm[22];
921:   if (npos) *npos = mat_mkl_pardiso->iparm[21];
922:   if (nzero) *nzero = F->rmap->N - (mat_mkl_pardiso->iparm[22] + mat_mkl_pardiso->iparm[21]);
923:   return(0);
924: }
925: #endif

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

932:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
933: #if defined(PETSC_USE_COMPLEX)
934:   F->ops->getinertia = NULL;
935: #else
936:   F->ops->getinertia = MatGetInertia_MKL_PARDISO;
937: #endif
938:   return(0);
939: }

941: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
942: {
943:   PetscErrorCode    ierr;
944:   PetscBool         iascii;
945:   PetscViewerFormat format;
946:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
947:   PetscInt          i;

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

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


973: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
974: {
975:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)A->data;

978:   info->block_size        = 1.0;
979:   info->nz_used           = mat_mkl_pardiso->iparm[17];
980:   info->nz_allocated      = mat_mkl_pardiso->iparm[17];
981:   info->nz_unneeded       = 0.0;
982:   info->assemblies        = 0.0;
983:   info->mallocs           = 0.0;
984:   info->memory            = 0.0;
985:   info->fill_ratio_given  = 0;
986:   info->fill_ratio_needed = 0;
987:   info->factor_mallocs    = 0;
988:   return(0);
989: }

991: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
992: {
993:   PetscInt        backup,bs;
994:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;

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

1022: /*@
1023:   MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters

1025:    Logically Collective on Mat

1027:    Input Parameters:
1028: +  F - the factored matrix obtained by calling MatGetFactor()
1029: .  icntl - index of Mkl_Pardiso parameter
1030: -  ival - value of Mkl_Pardiso parameter

1032:   Options Database:
1033: .   -mat_mkl_pardiso_<icntl> <ival>

1035:    Level: beginner

1037:    References:
1038: .      Mkl_Pardiso Users' Guide

1040: .seealso: MatGetFactor()
1041: @*/
1042: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
1043: {

1047:   PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
1048:   return(0);
1049: }

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

1055:   Works with MATSEQAIJ matrices

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

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

1085:   Level: beginner

1087:   For more information please check  mkl_pardiso manual

1089: .seealso: PCFactorSetMatSolverType(), MatSolverType

1091: M*/
1092: static PetscErrorCode MatFactorGetSolverType_mkl_pardiso(Mat A, MatSolverType *type)
1093: {
1095:   *type = MATSOLVERMKL_PARDISO;
1096:   return(0);
1097: }

1099: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
1100: {
1101:   Mat             B;
1102:   PetscErrorCode  ierr;
1103:   Mat_MKL_PARDISO *mat_mkl_pardiso;
1104:   PetscBool       isSeqAIJ,isSeqBAIJ,isSeqSBAIJ;

1107:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1108:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1109:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1110:   MatCreate(PetscObjectComm((PetscObject)A),&B);
1111:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1112:   PetscStrallocpy("mkl_pardiso",&((PetscObject)B)->type_name);
1113:   MatSetUp(B);

1115:   PetscNewLog(B,&mat_mkl_pardiso);
1116:   B->data = mat_mkl_pardiso;

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

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

1155:   PetscFree(B->solvertype);
1156:   PetscStrallocpy(MATSOLVERMKL_PARDISO,&B->solvertype);

1158:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_pardiso);
1159:   PetscObjectComposeFunction((PetscObject)B,"MatFactorSetSchurIS_C",MatFactorSetSchurIS_MKL_PARDISO);
1160:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);

1162:   *F = B;
1163:   return(0);
1164: }

1166: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void)
1167: {

1171:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1172:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1173:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQBAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1174:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1175:   return(0);
1176: }