Actual source code: mkl_pardiso.c

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

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

     10: PETSC_EXTERN void PetscSetMKL_PARDISOThreads(int);

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

     28: #define IPARM_SIZE 64

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    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:   VecRestoreArrayWrite(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:     MatDenseGetArrayWrite(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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    778:   PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
    779:   if (flg) mat_mkl_pardiso->iparm[59] = icntl;
    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 rank-0 */

    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:   PetscCalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
    844:   mat_mkl_pardiso->schur_size = 0;
    845:   return(0);
    846: }

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

    854:   mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
    855:   PetscSetMKL_PARDISOFromOptions(F,A);
    856:   /* throw away any previously computed structure */
    857:   if (mat_mkl_pardiso->freeaij) {
    858:     PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
    859:     if (mat_mkl_pardiso->iparm[34] == 1) {
    860:       PetscFree(mat_mkl_pardiso->a);
    861:     }
    862:   }
    863:   (*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);
    864:   if (mat_mkl_pardiso->iparm[34] == 1) mat_mkl_pardiso->n = A->rmap->N;
    865:   else mat_mkl_pardiso->n = A->rmap->N/A->rmap->bs;

    867:   mat_mkl_pardiso->phase = JOB_ANALYSIS;

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

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

    890:   mat_mkl_pardiso->CleanUp = PETSC_TRUE;

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

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

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

    906:   MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
    907:   return(0);
    908: }

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

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

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

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

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

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

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

    969: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
    970: {
    971:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)A->data;

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

    987: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
    988: {
    989:   PetscInt        backup,bs;
    990:   Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;

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

    1018: /*@
    1019:   MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters

    1021:    Logically Collective on Mat

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

    1028:   Options Database:
    1029: .   -mat_mkl_pardiso_<icntl> <ival>

    1031:    Level: beginner

    1033:    References:
    1034: .      Mkl_Pardiso Users' Guide

    1036: .seealso: MatGetFactor()
    1037: @*/
    1038: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
    1039: {

    1043:   PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
    1044:   return(0);
    1045: }

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

    1051:   Works with MATSEQAIJ matrices

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

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

    1081:   Level: beginner

    1083:   For more information please check  mkl_pardiso manual

    1085: .seealso: PCFactorSetMatSolverType(), MatSolverType

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

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

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

    1111:   PetscNewLog(B,&mat_mkl_pardiso);
    1112:   B->data = mat_mkl_pardiso;

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

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

    1151:   PetscFree(B->solvertype);
    1152:   PetscStrallocpy(MATSOLVERMKL_PARDISO,&B->solvertype);

    1154:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_pardiso);
    1155:   PetscObjectComposeFunction((PetscObject)B,"MatFactorSetSchurIS_C",MatFactorSetSchurIS_MKL_PARDISO);
    1156:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);

    1158:   *F = B;
    1159:   return(0);
    1160: }

    1162: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void)
    1163: {

    1167:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
    1168:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
    1169:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQBAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
    1170:   MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
    1171:   return(0);
    1172: }