Actual source code: mkl_pardiso.c

petsc-master 2019-06-22
Report Typos and Errors
  1: #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64)
  2: #define MKL_ILP64
  3: #endif

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

  9: #include <stdio.h>
 10: #include <stdlib.h>
 11: #include <math.h>
 12: #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64)
 13: #define MKL_ILP64
 14: #endif
 15: #include <mkl_pardiso.h>

 17: PETSC_EXTERN void PetscSetMKL_PARDISOThreads(int);

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

 35: #define IPARM_SIZE 64

 37: #if defined(PETSC_USE_64BIT_INDICES)
 38:  #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64)
 39:   #define INT_TYPE long long int
 40:   #define MKL_PARDISO pardiso
 41:   #define MKL_PARDISO_INIT pardisoinit
 42:  #else
 43:   /* this is the case where the MKL BLAS/LAPACK are 32 bit integers but the 64 bit integer version of 
 44:      of Pardiso code is used; hence the need for the 64 below*/
 45:   #define INT_TYPE long long int
 46:   #define MKL_PARDISO pardiso_64
 47:   #define MKL_PARDISO_INIT pardiso_64init
 48: void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm [])
 49: {
 50:   int iparm_copy[IPARM_SIZE], mtype_copy, i;
 51: 
 52:   mtype_copy = *mtype;
 53:   pardisoinit(pt, &mtype_copy, iparm_copy);
 54:   for(i = 0; i < IPARM_SIZE; i++){
 55:     iparm[i] = iparm_copy[i];
 56:   }
 57: }
 58:  #endif
 59: #else
 60:  #define INT_TYPE int
 61:  #define MKL_PARDISO pardiso
 62:  #define MKL_PARDISO_INIT pardisoinit
 63: #endif


 66: /*
 67:  *  Internal data structure.
 68:  *  For more information check mkl_pardiso manual.
 69:  */
 70: typedef struct {

 72:   /* Configuration vector*/
 73:   INT_TYPE     iparm[IPARM_SIZE];

 75:   /*
 76:    * Internal mkl_pardiso memory location.
 77:    * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
 78:    */
 79:   void         *pt[IPARM_SIZE];

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

 84:   /* Matrix structure*/
 85:   void         *a;
 86:   INT_TYPE     *ia, *ja;

 88:   /* Number of non-zero elements*/
 89:   INT_TYPE     nz;

 91:   /* Row permutaton vector*/
 92:   INT_TYPE     *perm;

 94:   /* Define if matrix preserves sparse structure.*/
 95:   MatStructure matstruc;

 97:   PetscBool    needsym;
 98:   PetscBool    freeaij;

100:   /* Schur complement */
101:   PetscScalar  *schur;
102:   PetscInt     schur_size;
103:   PetscInt     *schur_idxs;
104:   PetscScalar  *schur_work;
105:   PetscBLASInt schur_work_size;
106:   PetscBool    solve_interior;

108:   /* True if mkl_pardiso function have been used.*/
109:   PetscBool CleanUp;

111:   /* Conversion to a format suitable for MKL */
112:   PetscErrorCode (*Convert)(Mat, PetscBool, MatReuse, PetscBool*, INT_TYPE*, INT_TYPE**, INT_TYPE**, PetscScalar**);
113: } Mat_MKL_PARDISO;

115: PetscErrorCode MatMKLPardiso_Convert_seqsbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
116: {
117:   Mat_SeqSBAIJ   *aa = (Mat_SeqSBAIJ*)A->data;
118:   PetscInt       bs  = A->rmap->bs;

121:   if (!sym) {
122:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
123:   }
124:   if (bs == 1) { /* already in the correct format */
125:     *v    = aa->a;
126:     /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
127:     *r    = (INT_TYPE*)aa->i;
128:     *c    = (INT_TYPE*)aa->j;
129:     *nnz  = (INT_TYPE)aa->nz;
130:     *free = PETSC_FALSE;
131:   } else {
132:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Conversion from SeqSBAIJ to MKL Pardiso format still need to be implemented");
133:   }
134:   return(0);
135: }

137: PetscErrorCode MatMKLPardiso_Convert_seqbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
138: {
140:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Conversion from SeqBAIJ to MKL Pardiso format still need to be implemented");
141:   return(0);
142: }

144: PetscErrorCode MatMKLPardiso_Convert_seqaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
145: {
146:   Mat_SeqAIJ     *aa = (Mat_SeqAIJ*)A->data;

150:   if (!sym) { /* already in the correct format */
151:     *v    = aa->a;
152:     *r    = (INT_TYPE*)aa->i;
153:     *c    = (INT_TYPE*)aa->j;
154:     *nnz  = (INT_TYPE)aa->nz;
155:     *free = PETSC_FALSE;
156:     return(0);
157:   }
158:   /* need to get the triangular part */
159:   if (reuse == MAT_INITIAL_MATRIX) {
160:     PetscScalar *vals,*vv;
161:     PetscInt    *row,*col,*jj;
162:     PetscInt    m = A->rmap->n,nz,i;

164:     nz = 0;
165:     for (i=0; i<m; i++) {
166:       nz += aa->i[i+1] - aa->diag[i];
167:     }
168:     PetscMalloc2(m+1,&row,nz,&col);
169:     PetscMalloc1(nz,&vals);
170:     jj = col;
171:     vv = vals;

173:     row[0] = 0;
174:     for (i=0; i<m; i++) {
175:       PetscInt    *aj = aa->j + aa->diag[i];
176:       PetscScalar *av = aa->a + aa->diag[i];
177:       PetscInt    rl = aa->i[i+1] - aa->diag[i],j;
178:       for (j=0; j<rl; j++) {
179:         *jj = *aj; jj++; aj++;
180:         *vv = *av; vv++; av++;
181:       }
182:       row[i+1]    = row[i] + rl;
183:     }
184:     *v    = vals;
185:     *r    = (INT_TYPE*)row;
186:     *c    = (INT_TYPE*)col;
187:     *nnz  = (INT_TYPE)nz;
188:     *free = PETSC_TRUE;
189:   } else {
190:     PetscScalar *vv;
191:     PetscInt    m = A->rmap->n,i;

193:     vv = *v;
194:     for (i=0; i<m; i++) {
195:       PetscScalar *av = aa->a + aa->diag[i];
196:       PetscInt    rl = aa->i[i+1] - aa->diag[i],j;
197:       for (j=0; j<rl; j++) {
198:         *vv = *av; vv++; av++;
199:       }
200:     }
201:     *free = PETSC_TRUE;
202:   }
203:   return(0);
204: }


207: static PetscErrorCode MatMKLPardisoSolveSchur_Private(Mat F, PetscScalar *B, PetscScalar *X)
208: {
209:   Mat_MKL_PARDISO      *mpardiso =(Mat_MKL_PARDISO*)F->data;
210:   Mat                  S,Xmat,Bmat;
211:   MatFactorSchurStatus schurstatus;
212:   PetscErrorCode       ierr;

215:   MatFactorFactorizeSchurComplement(F);
216:   MatFactorGetSchurComplement(F,&S,&schurstatus);
217:   if (X == B && schurstatus == MAT_FACTOR_SCHUR_INVERTED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"X and B cannot point to the same address");
218:   MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,B,&Bmat);
219:   MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,X,&Xmat);
220:   if (X != B) { /* using MatMatSolve */
221:     MatCopy(Bmat,Xmat,SAME_NONZERO_PATTERN);
222:   }

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

228:   switch (schurstatus) {
229:   case MAT_FACTOR_SCHUR_FACTORED:
230:     if (!mpardiso->iparm[12-1]) {
231:       MatMatSolve(S,Bmat,Xmat);
232:     } else { /* transpose solve */
233:       MatMatSolveTranspose(S,Bmat,Xmat);
234:     }
235:     break;
236:   case MAT_FACTOR_SCHUR_INVERTED:
237:     if (!mpardiso->iparm[12-1]) {
238:       MatMatMult(S,Bmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Xmat);
239:     } else { /* transpose solve */
240:       MatTransposeMatMult(S,Bmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Xmat);
241:     }
242:     break;
243:   default:
244:     SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status);
245:     break;
246:   }
247:   MatFactorRestoreSchurComplement(F,&S,schurstatus);
248:   MatDestroy(&Bmat);
249:   MatDestroy(&Xmat);
250:   return(0);
251: }

253: PetscErrorCode MatFactorSetSchurIS_MKL_PARDISO(Mat F, IS is)
254: {
255:   Mat_MKL_PARDISO *mpardiso =(Mat_MKL_PARDISO*)F->data;
256:   const PetscInt  *idxs;
257:   PetscInt        size,i;
258:   PetscMPIInt     csize;
259:   PetscBool       sorted;
260:   PetscErrorCode  ierr;

263:   MPI_Comm_size(PetscObjectComm((PetscObject)F),&csize);
264:   if (csize > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MKL_PARDISO parallel Schur complements not yet supported from PETSc");
265:   ISSorted(is,&sorted);
266:   if (!sorted) {
267:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS for MKL_PARDISO Schur complements needs to be sorted");
268:   }
269:   ISGetLocalSize(is,&size);
270:   if (mpardiso->schur_size != size) {
271:     mpardiso->schur_size = size;
272:     PetscFree2(mpardiso->schur,mpardiso->schur_work);
273:     PetscFree(mpardiso->schur_idxs);
274:     PetscBLASIntCast(PetscMax(mpardiso->n,2*size),&mpardiso->schur_work_size);
275:     PetscMalloc2(size*size,&mpardiso->schur,mpardiso->schur_work_size,&mpardiso->schur_work);
276:     PetscMalloc1(size,&mpardiso->schur_idxs);
277:   }
278:   MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->schur_size,mpardiso->schur,&F->schur);
279:   if (mpardiso->mtype == 2) {
280:     MatSetOption(F->schur,MAT_SPD,PETSC_TRUE);
281:   }

283:   PetscArrayzero(mpardiso->perm,mpardiso->n);
284:   ISGetIndices(is,&idxs);
285:   PetscArraycpy(mpardiso->schur_idxs,idxs,size);
286:   for (i=0;i<size;i++) mpardiso->perm[idxs[i]] = 1;
287:   ISRestoreIndices(is,&idxs);
288:   if (size) { /* turn on Schur switch if the set of indices is not empty */
289:     mpardiso->iparm[36-1] = 2;
290:   }
291:   return(0);
292: }

294: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
295: {
296:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
297:   PetscErrorCode  ierr;

300:   if (mat_mkl_pardiso->CleanUp) {
301:     mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

303:     MKL_PARDISO (mat_mkl_pardiso->pt,
304:       &mat_mkl_pardiso->maxfct,
305:       &mat_mkl_pardiso->mnum,
306:       &mat_mkl_pardiso->mtype,
307:       &mat_mkl_pardiso->phase,
308:       &mat_mkl_pardiso->n,
309:       NULL,
310:       NULL,
311:       NULL,
312:       NULL,
313:       &mat_mkl_pardiso->nrhs,
314:       mat_mkl_pardiso->iparm,
315:       &mat_mkl_pardiso->msglvl,
316:       NULL,
317:       NULL,
318:       &mat_mkl_pardiso->err);
319:   }
320:   PetscFree(mat_mkl_pardiso->perm);
321:   PetscFree2(mat_mkl_pardiso->schur,mat_mkl_pardiso->schur_work);
322:   PetscFree(mat_mkl_pardiso->schur_idxs);
323:   if (mat_mkl_pardiso->freeaij) {
324:     PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
325:     PetscFree(mat_mkl_pardiso->a);
326:   }
327:   PetscFree(A->data);

329:   /* clear composed functions */
330:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
331:   PetscObjectComposeFunction((PetscObject)A,"MatFactorSetSchurIS_C",NULL);
332:   PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);
333:   return(0);
334: }

336: static PetscErrorCode MatMKLPardisoScatterSchur_Private(Mat_MKL_PARDISO *mpardiso, PetscScalar *whole, PetscScalar *schur, PetscBool reduce)
337: {
339:   if (reduce) { /* data given for the whole matrix */
340:     PetscInt i,m=0,p=0;
341:     for (i=0;i<mpardiso->nrhs;i++) {
342:       PetscInt j;
343:       for (j=0;j<mpardiso->schur_size;j++) {
344:         schur[p+j] = whole[m+mpardiso->schur_idxs[j]];
345:       }
346:       m += mpardiso->n;
347:       p += mpardiso->schur_size;
348:     }
349:   } else { /* from Schur to whole */
350:     PetscInt i,m=0,p=0;
351:     for (i=0;i<mpardiso->nrhs;i++) {
352:       PetscInt j;
353:       for (j=0;j<mpardiso->schur_size;j++) {
354:         whole[m+mpardiso->schur_idxs[j]] = schur[p+j];
355:       }
356:       m += mpardiso->n;
357:       p += mpardiso->schur_size;
358:     }
359:   }
360:   return(0);
361: }

363: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x)
364: {
365:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
366:   PetscErrorCode    ierr;
367:   PetscScalar       *xarray;
368:   const PetscScalar *barray;

371:   mat_mkl_pardiso->nrhs = 1;
372:   VecGetArray(x,&xarray);
373:   VecGetArrayRead(b,&barray);

375:   if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
376:   else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;

378:   if (barray == xarray) { /* if the two vectors share the same memory */
379:     PetscScalar *work;
380:     if (!mat_mkl_pardiso->schur_work) {
381:       PetscMalloc1(mat_mkl_pardiso->n,&work);
382:     } else {
383:       work = mat_mkl_pardiso->schur_work;
384:     }
385:     mat_mkl_pardiso->iparm[6-1] = 1;
386:     MKL_PARDISO (mat_mkl_pardiso->pt,
387:       &mat_mkl_pardiso->maxfct,
388:       &mat_mkl_pardiso->mnum,
389:       &mat_mkl_pardiso->mtype,
390:       &mat_mkl_pardiso->phase,
391:       &mat_mkl_pardiso->n,
392:       mat_mkl_pardiso->a,
393:       mat_mkl_pardiso->ia,
394:       mat_mkl_pardiso->ja,
395:       NULL,
396:       &mat_mkl_pardiso->nrhs,
397:       mat_mkl_pardiso->iparm,
398:       &mat_mkl_pardiso->msglvl,
399:       (void*)xarray,
400:       (void*)work,
401:       &mat_mkl_pardiso->err);
402:     if (!mat_mkl_pardiso->schur_work) {
403:       PetscFree(work);
404:     }
405:   } else {
406:     mat_mkl_pardiso->iparm[6-1] = 0;
407:     MKL_PARDISO (mat_mkl_pardiso->pt,
408:       &mat_mkl_pardiso->maxfct,
409:       &mat_mkl_pardiso->mnum,
410:       &mat_mkl_pardiso->mtype,
411:       &mat_mkl_pardiso->phase,
412:       &mat_mkl_pardiso->n,
413:       mat_mkl_pardiso->a,
414:       mat_mkl_pardiso->ia,
415:       mat_mkl_pardiso->ja,
416:       mat_mkl_pardiso->perm,
417:       &mat_mkl_pardiso->nrhs,
418:       mat_mkl_pardiso->iparm,
419:       &mat_mkl_pardiso->msglvl,
420:       (void*)barray,
421:       (void*)xarray,
422:       &mat_mkl_pardiso->err);
423:   }
424:   VecRestoreArrayRead(b,&barray);

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

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

431:     /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
432:     if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) {
433:       shift = 0;
434:     }

436:     if (!mat_mkl_pardiso->solve_interior) {
437:       /* solve Schur complement */
438:       MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
439:       MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
440:       MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
441:     } else { /* if we are solving for the interior problem, any value in barray[schur] forward-substituted to xarray[schur] will be neglected */
442:       PetscInt i;
443:       for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
444:         xarray[mat_mkl_pardiso->schur_idxs[i]] = 0.;
445:       }
446:     }

448:     /* expansion phase */
449:     mat_mkl_pardiso->iparm[6-1] = 1;
450:     mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
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*)xarray,
465:       (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
466:       &mat_mkl_pardiso->err);

468:     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:     mat_mkl_pardiso->iparm[6-1] = 0;
470:   }
471:   VecRestoreArray(x,&xarray);
472:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
473:   return(0);
474: }

476: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x)
477: {
478:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
479:   PetscInt        oiparm12;
480:   PetscErrorCode  ierr;

483:   oiparm12 = mat_mkl_pardiso->iparm[12 - 1];
484:   mat_mkl_pardiso->iparm[12 - 1] = 2;
485:   MatSolve_MKL_PARDISO(A,b,x);
486:   mat_mkl_pardiso->iparm[12 - 1] = oiparm12;
487:   return(0);
488: }

490: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X)
491: {
492:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->data;
493:   PetscErrorCode    ierr;
494:   const PetscScalar *barray;
495:   PetscScalar       *xarray;
496:   PetscBool         flg;

499:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
500:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
501:   PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
502:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");

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

506:   if (mat_mkl_pardiso->nrhs > 0) {
507:     MatDenseGetArrayRead(B,&barray);
508:     MatDenseGetArray(X,&xarray);

510:     if (barray == xarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"B and X cannot share the same memory location");
511:     if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
512:     else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
513:     mat_mkl_pardiso->iparm[6-1] = 0;

515:     MKL_PARDISO (mat_mkl_pardiso->pt,
516:       &mat_mkl_pardiso->maxfct,
517:       &mat_mkl_pardiso->mnum,
518:       &mat_mkl_pardiso->mtype,
519:       &mat_mkl_pardiso->phase,
520:       &mat_mkl_pardiso->n,
521:       mat_mkl_pardiso->a,
522:       mat_mkl_pardiso->ia,
523:       mat_mkl_pardiso->ja,
524:       mat_mkl_pardiso->perm,
525:       &mat_mkl_pardiso->nrhs,
526:       mat_mkl_pardiso->iparm,
527:       &mat_mkl_pardiso->msglvl,
528:       (void*)barray,
529:       (void*)xarray,
530:       &mat_mkl_pardiso->err);
531:     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);

533:     MatDenseRestoreArrayRead(B,&barray);
534:     if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
535:       PetscScalar *o_schur_work = NULL;
536:       PetscInt    shift = mat_mkl_pardiso->schur_size*mat_mkl_pardiso->nrhs,scale;
537:       PetscInt    mem = mat_mkl_pardiso->n*mat_mkl_pardiso->nrhs;

539:       /* allocate extra memory if it is needed */
540:       scale = 1;
541:       if (A->schur_status == MAT_FACTOR_SCHUR_INVERTED) scale = 2;

543:       mem *= scale;
544:       if (mem > mat_mkl_pardiso->schur_work_size) {
545:         o_schur_work = mat_mkl_pardiso->schur_work;
546:         PetscMalloc1(mem,&mat_mkl_pardiso->schur_work);
547:       }

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

552:       /* solve Schur complement */
553:       if (!mat_mkl_pardiso->solve_interior) {
554:         MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
555:         MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
556:         MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
557:       } else { /* if we are solving for the interior problem, any value in barray[schur,n] forward-substituted to xarray[schur,n] will be neglected */
558:         PetscInt i,n,m=0;
559:         for (n=0;n<mat_mkl_pardiso->nrhs;n++) {
560:           for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
561:             xarray[mat_mkl_pardiso->schur_idxs[i]+m] = 0.;
562:           }
563:           m += mat_mkl_pardiso->n;
564:         }
565:       }

567:       /* expansion phase */
568:       mat_mkl_pardiso->iparm[6-1] = 1;
569:       mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
570:       MKL_PARDISO (mat_mkl_pardiso->pt,
571:         &mat_mkl_pardiso->maxfct,
572:         &mat_mkl_pardiso->mnum,
573:         &mat_mkl_pardiso->mtype,
574:         &mat_mkl_pardiso->phase,
575:         &mat_mkl_pardiso->n,
576:         mat_mkl_pardiso->a,
577:         mat_mkl_pardiso->ia,
578:         mat_mkl_pardiso->ja,
579:         mat_mkl_pardiso->perm,
580:         &mat_mkl_pardiso->nrhs,
581:         mat_mkl_pardiso->iparm,
582:         &mat_mkl_pardiso->msglvl,
583:         (void*)xarray,
584:         (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
585:         &mat_mkl_pardiso->err);
586:       if (o_schur_work) { /* restore original schur_work (minimal size) */
587:         PetscFree(mat_mkl_pardiso->schur_work);
588:         mat_mkl_pardiso->schur_work = o_schur_work;
589:       }
590:       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);
591:       mat_mkl_pardiso->iparm[6-1] = 0;
592:     }
593:   }
594:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;
595:   return(0);
596: }

598: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info)
599: {
600:   Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->data;
601:   PetscErrorCode  ierr;

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

607:   mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
608:   MKL_PARDISO (mat_mkl_pardiso->pt,
609:     &mat_mkl_pardiso->maxfct,
610:     &mat_mkl_pardiso->mnum,
611:     &mat_mkl_pardiso->mtype,
612:     &mat_mkl_pardiso->phase,
613:     &mat_mkl_pardiso->n,
614:     mat_mkl_pardiso->a,
615:     mat_mkl_pardiso->ia,
616:     mat_mkl_pardiso->ja,
617:     mat_mkl_pardiso->perm,
618:     &mat_mkl_pardiso->nrhs,
619:     mat_mkl_pardiso->iparm,
620:     &mat_mkl_pardiso->msglvl,
621:     NULL,
622:     (void*)mat_mkl_pardiso->schur,
623:     &mat_mkl_pardiso->err);
624:   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);

626:   if (F->schur) { /* schur output from pardiso is in row major format */
627:     MatFactorRestoreSchurComplement(F,NULL,MAT_FACTOR_SCHUR_UNFACTORED);
628:     MatTranspose(F->schur,MAT_INPLACE_MATRIX,&F->schur);
629:   }
630:   mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
631:   mat_mkl_pardiso->CleanUp  = PETSC_TRUE;
632:   return(0);
633: }

635: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A)
636: {
637:   Mat_MKL_PARDISO     *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
638:   PetscErrorCode      ierr;
639:   PetscInt            icntl,threads=1;
640:   PetscBool           flg;

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

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

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

651:   PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);
652:   if (flg) mat_mkl_pardiso->mnum = icntl;
653: 
654:   PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);
655:   if (flg) mat_mkl_pardiso->msglvl = icntl;

657:   PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);
658:   if(flg){
659:     void *pt[IPARM_SIZE];
660:     mat_mkl_pardiso->mtype = icntl;
661:     MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
662: #if defined(PETSC_USE_REAL_SINGLE)
663:     mat_mkl_pardiso->iparm[27] = 1;
664: #else
665:     mat_mkl_pardiso->iparm[27] = 0;
666: #endif
667:     mat_mkl_pardiso->iparm[34] = 1; /* use 0-based indexing */
668:   }
669:   PetscOptionsInt("-mat_mkl_pardiso_1","Use default values","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

723:     PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
724:     if (flg) mat_mkl_pardiso->iparm[59] = icntl;
725:   }
726:   PetscOptionsEnd();
727:   return(0);
728: }

730: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
731: {
733:   PetscInt       i;

736:   for ( i = 0; i < IPARM_SIZE; i++ ){
737:     mat_mkl_pardiso->iparm[i] = 0;
738:   }
739:   for ( i = 0; i < IPARM_SIZE; i++ ){
740:     mat_mkl_pardiso->pt[i] = 0;
741:   }
742:   /* Default options for both sym and unsym */
743:   mat_mkl_pardiso->iparm[ 0] =  1; /* Solver default parameters overriden with provided by iparm */
744:   mat_mkl_pardiso->iparm[ 1] =  2; /* Metis reordering */
745:   mat_mkl_pardiso->iparm[ 5] =  0; /* Write solution into x */
746:   mat_mkl_pardiso->iparm[ 7] =  0; /* Max number of iterative refinement steps */
747:   mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
748:   mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
749: #if 0
750:   mat_mkl_pardiso->iparm[23] =  1; /* Parallel factorization control*/
751: #endif
752:   mat_mkl_pardiso->iparm[34] =  1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
753:   mat_mkl_pardiso->iparm[39] =  0; /* Input: matrix/rhs/solution stored on master */
754: 
755:   mat_mkl_pardiso->CleanUp   = PETSC_FALSE;
756:   mat_mkl_pardiso->maxfct    = 1; /* Maximum number of numerical factorizations. */
757:   mat_mkl_pardiso->mnum      = 1; /* Which factorization to use. */
758:   mat_mkl_pardiso->msglvl    = 0; /* 0: do not print 1: Print statistical information in file */
759:   mat_mkl_pardiso->phase     = -1;
760:   mat_mkl_pardiso->err       = 0;
761: 
762:   mat_mkl_pardiso->n         = A->rmap->N;
763:   mat_mkl_pardiso->nrhs      = 1;
764:   mat_mkl_pardiso->err       = 0;
765:   mat_mkl_pardiso->phase     = -1;
766: 
767:   if(ftype == MAT_FACTOR_LU){
768:     mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
769:     mat_mkl_pardiso->iparm[10] =  1; /* Use nonsymmetric permutation and scaling MPS */
770:     mat_mkl_pardiso->iparm[12] =  1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */

772:   } else {
773:     mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
774:     mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
775:     mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
776: /*    mat_mkl_pardiso->iparm[20] =  1; */ /* Apply 1x1 and 2x2 Bunch-Kaufman pivoting during the factorization process */
777: #if defined(PETSC_USE_DEBUG)
778:     mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
779: #endif
780:   }
781:   PetscMalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
782:   for(i = 0; i < A->rmap->N; i++){
783:     mat_mkl_pardiso->perm[i] = 0;
784:   }
785:   mat_mkl_pardiso->schur_size = 0;
786:   return(0);
787: }

789: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
790: {
791:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
792:   PetscErrorCode  ierr;

795:   mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
796:   PetscSetMKL_PARDISOFromOptions(F,A);

798:   /* throw away any previously computed structure */
799:   if (mat_mkl_pardiso->freeaij) {
800:     PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
801:     PetscFree(mat_mkl_pardiso->a);
802:   }
803:   (*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);
804:   mat_mkl_pardiso->n = A->rmap->N;

806:   mat_mkl_pardiso->phase = JOB_ANALYSIS;

808:   MKL_PARDISO (mat_mkl_pardiso->pt,
809:     &mat_mkl_pardiso->maxfct,
810:     &mat_mkl_pardiso->mnum,
811:     &mat_mkl_pardiso->mtype,
812:     &mat_mkl_pardiso->phase,
813:     &mat_mkl_pardiso->n,
814:     mat_mkl_pardiso->a,
815:     mat_mkl_pardiso->ia,
816:     mat_mkl_pardiso->ja,
817:     mat_mkl_pardiso->perm,
818:     &mat_mkl_pardiso->nrhs,
819:     mat_mkl_pardiso->iparm,
820:     &mat_mkl_pardiso->msglvl,
821:     NULL,
822:     NULL,
823:     &mat_mkl_pardiso->err);
824:   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);

826:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;

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

831:   F->ops->solve           = MatSolve_MKL_PARDISO;
832:   F->ops->solvetranspose  = MatSolveTranspose_MKL_PARDISO;
833:   F->ops->matsolve        = MatMatSolve_MKL_PARDISO;
834:   return(0);
835: }

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

842:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
843:   return(0);
844: }

846: #if !defined(PETSC_USE_COMPLEX)
847: PetscErrorCode MatGetInertia_MKL_PARDISO(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
848: {
849:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)F->data;

852:   if (nneg) *nneg = mat_mkl_pardiso->iparm[22];
853:   if (npos) *npos = mat_mkl_pardiso->iparm[21];
854:   if (nzero) *nzero = F->rmap->N -(mat_mkl_pardiso->iparm[22] + mat_mkl_pardiso->iparm[21]);
855:   return(0);
856: }
857: #endif

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

864:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
865: #if defined(PETSC_USE_COMPLEX)
866:   F->ops->getinertia = NULL;
867: #else
868:   F->ops->getinertia = MatGetInertia_MKL_PARDISO;
869: #endif
870:   return(0);
871: }

873: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
874: {
875:   PetscErrorCode    ierr;
876:   PetscBool         iascii;
877:   PetscViewerFormat format;
878:   Mat_MKL_PARDISO   *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
879:   PetscInt          i;

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

884:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
885:   if (iascii) {
886:     PetscViewerGetFormat(viewer,&format);
887:     if (format == PETSC_VIEWER_ASCII_INFO) {
888:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
889:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase:             %d \n",mat_mkl_pardiso->phase);
890:       for(i = 1; i <= 64; i++){
891:         PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]:     %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
892:       }
893:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct:     %d \n", mat_mkl_pardiso->maxfct);
894:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum:     %d \n", mat_mkl_pardiso->mnum);
895:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype:     %d \n", mat_mkl_pardiso->mtype);
896:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n:     %d \n", mat_mkl_pardiso->n);
897:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs:     %d \n", mat_mkl_pardiso->nrhs);
898:       PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl:     %d \n", mat_mkl_pardiso->msglvl);
899:     }
900:   }
901:   return(0);
902: }


905: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
906: {
907:   Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->data;

910:   info->block_size        = 1.0;
911:   info->nz_used           = mat_mkl_pardiso->nz;
912:   info->nz_allocated      = mat_mkl_pardiso->nz;
913:   info->nz_unneeded       = 0.0;
914:   info->assemblies        = 0.0;
915:   info->mallocs           = 0.0;
916:   info->memory            = 0.0;
917:   info->fill_ratio_given  = 0;
918:   info->fill_ratio_needed = 0;
919:   info->factor_mallocs    = 0;
920:   return(0);
921: }

923: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
924: {
925:   Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->data;

928:   if(icntl <= 64){
929:     mat_mkl_pardiso->iparm[icntl - 1] = ival;
930:   } else {
931:     if(icntl == 65) PetscSetMKL_PARDISOThreads(ival);
932:     else if(icntl == 66) mat_mkl_pardiso->maxfct = ival;
933:     else if(icntl == 67) mat_mkl_pardiso->mnum = ival;
934:     else if(icntl == 68) mat_mkl_pardiso->msglvl = ival;
935:     else if(icntl == 69){
936:       void *pt[IPARM_SIZE];
937:       mat_mkl_pardiso->mtype = ival;
938:       MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
939: #if defined(PETSC_USE_REAL_SINGLE)
940:       mat_mkl_pardiso->iparm[27] = 1;
941: #else
942:       mat_mkl_pardiso->iparm[27] = 0;
943: #endif
944:       mat_mkl_pardiso->iparm[34] = 1;
945:     } else if(icntl==70) mat_mkl_pardiso->solve_interior = (PetscBool)!!ival;
946:   }
947:   return(0);
948: }

950: /*@
951:   MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters

953:    Logically Collective on Mat

955:    Input Parameters:
956: +  F - the factored matrix obtained by calling MatGetFactor()
957: .  icntl - index of Mkl_Pardiso parameter
958: -  ival - value of Mkl_Pardiso parameter

960:   Options Database:
961: .   -mat_mkl_pardiso_<icntl> <ival>

963:    Level: beginner

965:    References:
966: .      Mkl_Pardiso Users' Guide

968: .seealso: MatGetFactor()
969: @*/
970: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
971: {

975:   PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
976:   return(0);
977: }

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

983:   Works with MATSEQAIJ matrices

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

987:   Options Database Keys:
988: + -mat_mkl_pardiso_65 - Number of threads to use within MKL_PARDISO
989: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
990: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
991: . -mat_mkl_pardiso_68 - Message level information
992: . -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
993: . -mat_mkl_pardiso_1  - Use default values
994: . -mat_mkl_pardiso_2  - Fill-in reducing ordering for the input matrix
995: . -mat_mkl_pardiso_4  - Preconditioned CGS/CG
996: . -mat_mkl_pardiso_5  - User permutation
997: . -mat_mkl_pardiso_6  - Write solution on x
998: . -mat_mkl_pardiso_8  - Iterative refinement step
999: . -mat_mkl_pardiso_10 - Pivoting perturbation
1000: . -mat_mkl_pardiso_11 - Scaling vectors
1001: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
1002: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
1003: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
1004: . -mat_mkl_pardiso_19 - Report number of floating point operations
1005: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
1006: . -mat_mkl_pardiso_24 - Parallel factorization control
1007: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
1008: . -mat_mkl_pardiso_27 - Matrix checker
1009: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
1010: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
1011: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode

1013:   Level: beginner

1015:   For more information please check  mkl_pardiso manual

1017: .seealso: PCFactorSetMatSolverType(), MatSolverType

1019: M*/
1020: static PetscErrorCode MatFactorGetSolverType_mkl_pardiso(Mat A, MatSolverType *type)
1021: {
1023:   *type = MATSOLVERMKL_PARDISO;
1024:   return(0);
1025: }

1027: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
1028: {
1029:   Mat             B;
1030:   PetscErrorCode  ierr;
1031:   Mat_MKL_PARDISO *mat_mkl_pardiso;
1032:   PetscBool       isSeqAIJ,isSeqBAIJ,isSeqSBAIJ;

1035:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1036:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1037:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1038:   MatCreate(PetscObjectComm((PetscObject)A),&B);
1039:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1040:   PetscStrallocpy("mkl_pardiso",&((PetscObject)B)->type_name);
1041:   MatSetUp(B);

1043:   PetscNewLog(B,&mat_mkl_pardiso);
1044:   B->data = mat_mkl_pardiso;

1046:   MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
1047:   if (ftype == MAT_FACTOR_LU) {
1048:     B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
1049:     B->factortype            = MAT_FACTOR_LU;
1050:     mat_mkl_pardiso->needsym = PETSC_FALSE;
1051:     if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1052:     else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1053:     else if (isSeqSBAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU factor with SEQSBAIJ format! Use MAT_FACTOR_CHOLESKY instead");
1054:     else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU with %s format",((PetscObject)A)->type_name);
1055: #if defined(PETSC_USE_COMPLEX)
1056:     mat_mkl_pardiso->mtype = 13;
1057: #else
1058:     mat_mkl_pardiso->mtype = 11;
1059: #endif
1060:   } else {
1061: #if defined(PETSC_USE_COMPLEX)
1062:     SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead",((PetscObject)A)->type_name);
1063: #endif
1064:     B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
1065:     B->factortype                  = MAT_FACTOR_CHOLESKY;
1066:     if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1067:     else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1068:     else if (isSeqSBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqsbaij;
1069:     else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with %s format",((PetscObject)A)->type_name);

1071:     mat_mkl_pardiso->needsym = PETSC_TRUE;
1072:     if (A->spd_set && A->spd) mat_mkl_pardiso->mtype = 2;
1073:     else                      mat_mkl_pardiso->mtype = -2;
1074:   }
1075:   B->ops->destroy          = MatDestroy_MKL_PARDISO;
1076:   B->ops->view             = MatView_MKL_PARDISO;
1077:   B->factortype            = ftype;
1078:   B->ops->getinfo          = MatGetInfo_MKL_PARDISO;
1079:   B->assembled             = PETSC_TRUE;

1081:   PetscFree(B->solvertype);
1082:   PetscStrallocpy(MATSOLVERMKL_PARDISO,&B->solvertype);

1084:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_pardiso);
1085:   PetscObjectComposeFunction((PetscObject)B,"MatFactorSetSchurIS_C",MatFactorSetSchurIS_MKL_PARDISO);
1086:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);

1088:   *F = B;
1089:   return(0);
1090: }

1092: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void)
1093: {

1097:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1098:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1099:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1100:   return(0);
1101: }