Actual source code: mkl_cpardiso.c

petsc-master 2019-07-18
Report Typos and Errors

  2: #include <../src/mat/impls/aij/seq/aij.h>                       /*I "petscmat.h" I*/
  3: #include <../src/mat/impls/aij/mpi/mpiaij.h>

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

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

 27: #define IPARM_SIZE 64
 28: #define INT_TYPE MKL_INT

 30: static const char *Err_MSG_CPardiso(int errNo){
 31:   switch (errNo) {
 32:     case -1:
 33:       return "input inconsistent"; break;
 34:     case -2:
 35:       return "not enough memory"; break;
 36:     case -3:
 37:       return "reordering problem"; break;
 38:     case -4:
 39:       return "zero pivot, numerical factorization or iterative refinement problem"; break;
 40:     case -5:
 41:       return "unclassified (internal) error"; break;
 42:     case -6:
 43:       return "preordering failed (matrix types 11, 13 only)"; break;
 44:     case -7:
 45:       return "diagonal matrix problem"; break;
 46:     case -8:
 47:       return "32-bit integer overflow problem"; break;
 48:     case -9:
 49:       return "not enough memory for OOC"; break;
 50:     case -10:
 51:       return "problems with opening OOC temporary files"; break;
 52:     case -11:
 53:       return "read/write problems with the OOC data file"; break;
 54:     default :
 55:       return "unknown error";
 56:   }
 57: }

 59: /*
 60:  *  Internal data structure.
 61:  *  For more information check mkl_cpardiso manual.
 62:  */

 64: typedef struct {

 66:   /* Configuration vector */
 67:   INT_TYPE     iparm[IPARM_SIZE];

 69:   /*
 70:    * Internal mkl_cpardiso memory location.
 71:    * After the first call to mkl_cpardiso do not modify pt, as that could cause a serious memory leak.
 72:    */
 73:   void         *pt[IPARM_SIZE];

 75:   MPI_Comm     comm_mkl_cpardiso;

 77:   /* Basic mkl_cpardiso info*/
 78:   INT_TYPE     phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;

 80:   /* Matrix structure */
 81:   PetscScalar  *a;

 83:   INT_TYPE     *ia, *ja;

 85:   /* Number of non-zero elements */
 86:   INT_TYPE     nz;

 88:   /* Row permutaton vector*/
 89:   INT_TYPE     *perm;

 91:   /* Define is matrix preserve sparce structure. */
 92:   MatStructure matstruc;

 94:   PetscErrorCode (*ConvertToTriples)(Mat, MatReuse, PetscInt*, PetscInt**, PetscInt**, PetscScalar**);

 96:   /* True if mkl_cpardiso function have been used. */
 97:   PetscBool CleanUp;
 98: } Mat_MKL_CPARDISO;

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

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

130: PetscErrorCode MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
131: {
132:   const PetscInt    *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
133:   PetscErrorCode    ierr;
134:   PetscInt          rstart,nz,i,j,countA,countB;
135:   PetscInt          *row,*col;
136:   const PetscScalar *av, *bv;
137:   PetscScalar       *val;
138:   Mat_MPIAIJ        *mat = (Mat_MPIAIJ*)A->data;
139:   Mat_SeqAIJ        *aa  = (Mat_SeqAIJ*)(mat->A)->data;
140:   Mat_SeqAIJ        *bb  = (Mat_SeqAIJ*)(mat->B)->data;
141:   PetscInt          colA_start,jB,jcol;

144:   ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart;
145:   av=aa->a; bv=bb->a;

147:   garray = mat->garray;

149:   if (reuse == MAT_INITIAL_MATRIX) {
150:     nz   = aa->nz + bb->nz;
151:     *nnz = nz;
152:     PetscMalloc((nz*(sizeof(PetscInt)+sizeof(PetscScalar)) + (m+1)*sizeof(PetscInt)), &row);
153:     col  = row + m + 1;
154:     val  = (PetscScalar*)(col + nz);
155:     *r = row; *c = col; *v = val;
156:     row[0] = 0;
157:   } else {
158:     row = *r; col = *c; val = *v;
159:   }

161:   nz = 0;
162:   for (i=0; i<m; i++) {
163:     row[i] = nz;
164:     countA     = ai[i+1] - ai[i];
165:     countB     = bi[i+1] - bi[i];
166:     ajj        = aj + ai[i]; /* ptr to the beginning of this row */
167:     bjj        = bj + bi[i];

169:     /* B part, smaller col index */
170:     colA_start = rstart + ajj[0]; /* the smallest global col index of A */
171:     jB         = 0;
172:     for (j=0; j<countB; j++) {
173:       jcol = garray[bjj[j]];
174:       if (jcol > colA_start) {
175:         jB = j;
176:         break;
177:       }
178:       col[nz]   = jcol;
179:       val[nz++] = *bv++;
180:       if (j==countB-1) jB = countB;
181:     }

183:     /* A part */
184:     for (j=0; j<countA; j++) {
185:       col[nz]   = rstart + ajj[j];
186:       val[nz++] = *av++;
187:     }

189:     /* B part, larger col index */
190:     for (j=jB; j<countB; j++) {
191:       col[nz]   = garray[bjj[j]];
192:       val[nz++] = *bv++;
193:     }
194:   }
195:   row[m] = nz;

197:   return(0);
198: }

200: /*
201:  * Free memory for Mat_MKL_CPARDISO structure and pointers to objects.
202:  */
203: PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A)
204: {
205:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
206:   PetscErrorCode   ierr;

209:   /* Terminate instance, deallocate memories */
210:   if (mat_mkl_cpardiso->CleanUp) {
211:     mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

213:     cluster_sparse_solver (
214:       mat_mkl_cpardiso->pt,
215:       &mat_mkl_cpardiso->maxfct,
216:       &mat_mkl_cpardiso->mnum,
217:       &mat_mkl_cpardiso->mtype,
218:       &mat_mkl_cpardiso->phase,
219:       &mat_mkl_cpardiso->n,
220:       NULL,
221:       NULL,
222:       NULL,
223:       mat_mkl_cpardiso->perm,
224:       &mat_mkl_cpardiso->nrhs,
225:       mat_mkl_cpardiso->iparm,
226:       &mat_mkl_cpardiso->msglvl,
227:       NULL,
228:       NULL,
229:       &mat_mkl_cpardiso->comm_mkl_cpardiso,
230:       (PetscInt*)&mat_mkl_cpardiso->err);
231:   }

233:   if (mat_mkl_cpardiso->ConvertToTriples == MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO) {
234:     PetscFree(mat_mkl_cpardiso->ia);
235:   }
236:   MPI_Comm_free(&(mat_mkl_cpardiso->comm_mkl_cpardiso));
237:   PetscFree(A->data);

239:   /* clear composed functions */
240:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
241:   PetscObjectComposeFunction((PetscObject)A,"MatMkl_CPardisoSetCntl_C",NULL);
242:   return(0);
243: }

245: /*
246:  * Computes Ax = b
247:  */
248: PetscErrorCode MatSolve_MKL_CPARDISO(Mat A,Vec b,Vec x)
249: {
250:   Mat_MKL_CPARDISO   *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
251:   PetscErrorCode    ierr;
252:   PetscScalar       *xarray;
253:   const PetscScalar *barray;

256:   mat_mkl_cpardiso->nrhs = 1;
257:   VecGetArray(x,&xarray);
258:   VecGetArrayRead(b,&barray);

260:   /* solve phase */
261:   /*-------------*/
262:   mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
263:   cluster_sparse_solver (
264:     mat_mkl_cpardiso->pt,
265:     &mat_mkl_cpardiso->maxfct,
266:     &mat_mkl_cpardiso->mnum,
267:     &mat_mkl_cpardiso->mtype,
268:     &mat_mkl_cpardiso->phase,
269:     &mat_mkl_cpardiso->n,
270:     mat_mkl_cpardiso->a,
271:     mat_mkl_cpardiso->ia,
272:     mat_mkl_cpardiso->ja,
273:     mat_mkl_cpardiso->perm,
274:     &mat_mkl_cpardiso->nrhs,
275:     mat_mkl_cpardiso->iparm,
276:     &mat_mkl_cpardiso->msglvl,
277:     (void*)barray,
278:     (void*)xarray,
279:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
280:     (PetscInt*)&mat_mkl_cpardiso->err);

282:   if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));

284:   VecRestoreArray(x,&xarray);
285:   VecRestoreArrayRead(b,&barray);
286:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
287:   return(0);
288: }

290: PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A,Vec b,Vec x)
291: {
292:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
293:   PetscErrorCode   ierr;

296: #if defined(PETSC_USE_COMPLEX)
297:   mat_mkl_cpardiso->iparm[12 - 1] = 1;
298: #else
299:   mat_mkl_cpardiso->iparm[12 - 1] = 2;
300: #endif
301:   MatSolve_MKL_CPARDISO(A,b,x);
302:   mat_mkl_cpardiso->iparm[12 - 1] = 0;
303:   return(0);
304: }

306: PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A,Mat B,Mat X)
307: {
308:   Mat_MKL_CPARDISO  *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
309:   PetscErrorCode    ierr;
310:   PetscScalar       *xarray;
311:   const PetscScalar *barray;

314:   MatGetSize(B,NULL,(PetscInt*)&mat_mkl_cpardiso->nrhs);

316:   if(mat_mkl_cpardiso->nrhs > 0){
317:     MatDenseGetArrayRead(B,&barray);
318:     MatDenseGetArray(X,&xarray);

320:     /* solve phase */
321:     /*-------------*/
322:     mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
323:     cluster_sparse_solver (
324:       mat_mkl_cpardiso->pt,
325:       &mat_mkl_cpardiso->maxfct,
326:       &mat_mkl_cpardiso->mnum,
327:       &mat_mkl_cpardiso->mtype,
328:       &mat_mkl_cpardiso->phase,
329:       &mat_mkl_cpardiso->n,
330:       mat_mkl_cpardiso->a,
331:       mat_mkl_cpardiso->ia,
332:       mat_mkl_cpardiso->ja,
333:       mat_mkl_cpardiso->perm,
334:       &mat_mkl_cpardiso->nrhs,
335:       mat_mkl_cpardiso->iparm,
336:       &mat_mkl_cpardiso->msglvl,
337:       (void*)barray,
338:       (void*)xarray,
339:       &mat_mkl_cpardiso->comm_mkl_cpardiso,
340:       (PetscInt*)&mat_mkl_cpardiso->err);
341:     if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));
342:     MatDenseRestoreArrayRead(B,&barray);
343:     MatDenseRestoreArray(X,&xarray);

345:   }
346:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
347:   return(0);

349: }

351: /*
352:  * LU Decomposition
353:  */
354: PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F,Mat A,const MatFactorInfo *info)
355: {
356:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(F)->data;
357:   PetscErrorCode   ierr;

360:   mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
361:   (*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_REUSE_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);

363:   mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION;
364:   cluster_sparse_solver (
365:     mat_mkl_cpardiso->pt,
366:     &mat_mkl_cpardiso->maxfct,
367:     &mat_mkl_cpardiso->mnum,
368:     &mat_mkl_cpardiso->mtype,
369:     &mat_mkl_cpardiso->phase,
370:     &mat_mkl_cpardiso->n,
371:     mat_mkl_cpardiso->a,
372:     mat_mkl_cpardiso->ia,
373:     mat_mkl_cpardiso->ja,
374:     mat_mkl_cpardiso->perm,
375:     &mat_mkl_cpardiso->nrhs,
376:     mat_mkl_cpardiso->iparm,
377:     &mat_mkl_cpardiso->msglvl,
378:     NULL,
379:     NULL,
380:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
381:     &mat_mkl_cpardiso->err);
382:   if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));

384:   mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
385:   mat_mkl_cpardiso->CleanUp  = PETSC_TRUE;
386:   return(0);
387: }

389: /* Sets mkl_cpardiso options from the options database */
390: PetscErrorCode PetscSetMKL_CPARDISOFromOptions(Mat F, Mat A)
391: {
392:   Mat_MKL_CPARDISO    *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
393:   PetscErrorCode      ierr;
394:   PetscInt            icntl,threads;
395:   PetscBool           flg;

398:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_CPARDISO Options","Mat");
399:   PetscOptionsInt("-mat_mkl_cpardiso_65","Number of threads to use","None",threads,&threads,&flg);
400:   if (flg) mkl_set_num_threads((int)threads);

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

405:   PetscOptionsInt("-mat_mkl_cpardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_cpardiso->mnum,&icntl,&flg);
406:   if (flg) mat_mkl_cpardiso->mnum = icntl;

408:   PetscOptionsInt("-mat_mkl_cpardiso_68","Message level information","None",mat_mkl_cpardiso->msglvl,&icntl,&flg);
409:   if (flg) mat_mkl_cpardiso->msglvl = icntl;

411:   PetscOptionsInt("-mat_mkl_cpardiso_69","Defines the matrix type","None",mat_mkl_cpardiso->mtype,&icntl,&flg);
412:   if(flg){
413:     mat_mkl_cpardiso->mtype = icntl;
414: #if defined(PETSC_USE_REAL_SINGLE)
415:     mat_mkl_cpardiso->iparm[27] = 1;
416: #else
417:     mat_mkl_cpardiso->iparm[27] = 0;
418: #endif
419:     mat_mkl_cpardiso->iparm[34] = 1;
420:   }
421:   PetscOptionsInt("-mat_mkl_cpardiso_1","Use default values","None",mat_mkl_cpardiso->iparm[0],&icntl,&flg);

423:   if(flg && icntl != 0){
424:     PetscOptionsInt("-mat_mkl_cpardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_cpardiso->iparm[1],&icntl,&flg);
425:     if (flg) mat_mkl_cpardiso->iparm[1] = icntl;

427:     PetscOptionsInt("-mat_mkl_cpardiso_4","Preconditioned CGS/CG","None",mat_mkl_cpardiso->iparm[3],&icntl,&flg);
428:     if (flg) mat_mkl_cpardiso->iparm[3] = icntl;

430:     PetscOptionsInt("-mat_mkl_cpardiso_5","User permutation","None",mat_mkl_cpardiso->iparm[4],&icntl,&flg);
431:     if (flg) mat_mkl_cpardiso->iparm[4] = icntl;

433:     PetscOptionsInt("-mat_mkl_cpardiso_6","Write solution on x","None",mat_mkl_cpardiso->iparm[5],&icntl,&flg);
434:     if (flg) mat_mkl_cpardiso->iparm[5] = icntl;

436:     PetscOptionsInt("-mat_mkl_cpardiso_8","Iterative refinement step","None",mat_mkl_cpardiso->iparm[7],&icntl,&flg);
437:     if (flg) mat_mkl_cpardiso->iparm[7] = icntl;

439:     PetscOptionsInt("-mat_mkl_cpardiso_10","Pivoting perturbation","None",mat_mkl_cpardiso->iparm[9],&icntl,&flg);
440:     if (flg) mat_mkl_cpardiso->iparm[9] = icntl;

442:     PetscOptionsInt("-mat_mkl_cpardiso_11","Scaling vectors","None",mat_mkl_cpardiso->iparm[10],&icntl,&flg);
443:     if (flg) mat_mkl_cpardiso->iparm[10] = icntl;

445:     PetscOptionsInt("-mat_mkl_cpardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_cpardiso->iparm[11],&icntl,&flg);
446:     if (flg) mat_mkl_cpardiso->iparm[11] = icntl;

448:     PetscOptionsInt("-mat_mkl_cpardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_cpardiso->iparm[12],&icntl,
449:       &flg);
450:     if (flg) mat_mkl_cpardiso->iparm[12] = icntl;

452:     PetscOptionsInt("-mat_mkl_cpardiso_18","Numbers of non-zero elements","None",mat_mkl_cpardiso->iparm[17],&icntl,
453:       &flg);
454:     if (flg) mat_mkl_cpardiso->iparm[17] = icntl;

456:     PetscOptionsInt("-mat_mkl_cpardiso_19","Report number of floating point operations","None",mat_mkl_cpardiso->iparm[18],&icntl,&flg);
457:     if (flg) mat_mkl_cpardiso->iparm[18] = icntl;

459:     PetscOptionsInt("-mat_mkl_cpardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_cpardiso->iparm[20],&icntl,&flg);
460:     if (flg) mat_mkl_cpardiso->iparm[20] = icntl;

462:     PetscOptionsInt("-mat_mkl_cpardiso_24","Parallel factorization control","None",mat_mkl_cpardiso->iparm[23],&icntl,&flg);
463:     if (flg) mat_mkl_cpardiso->iparm[23] = icntl;

465:     PetscOptionsInt("-mat_mkl_cpardiso_25","Parallel forward/backward solve control","None",mat_mkl_cpardiso->iparm[24],&icntl,&flg);
466:     if (flg) mat_mkl_cpardiso->iparm[24] = icntl;

468:     PetscOptionsInt("-mat_mkl_cpardiso_27","Matrix checker","None",mat_mkl_cpardiso->iparm[26],&icntl,&flg);
469:     if (flg) mat_mkl_cpardiso->iparm[26] = icntl;

471:     PetscOptionsInt("-mat_mkl_cpardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_cpardiso->iparm[30],&icntl,&flg);
472:     if (flg) mat_mkl_cpardiso->iparm[30] = icntl;

474:     PetscOptionsInt("-mat_mkl_cpardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_cpardiso->iparm[33],&icntl,&flg);
475:     if (flg) mat_mkl_cpardiso->iparm[33] = icntl;

477:     PetscOptionsInt("-mat_mkl_cpardiso_60","Intel MKL_CPARDISO mode","None",mat_mkl_cpardiso->iparm[59],&icntl,&flg);
478:     if (flg) mat_mkl_cpardiso->iparm[59] = icntl;
479:   }

481:   PetscOptionsEnd();
482:   return(0);
483: }

485: PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso)
486: {
487:   PetscErrorCode  ierr;
488:   PetscMPIInt     size;


492:   MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mat_mkl_cpardiso->comm_mkl_cpardiso));
493:   MPI_Comm_size(mat_mkl_cpardiso->comm_mkl_cpardiso, &size);

495:   mat_mkl_cpardiso->CleanUp = PETSC_FALSE;
496:   mat_mkl_cpardiso->maxfct = 1;
497:   mat_mkl_cpardiso->mnum = 1;
498:   mat_mkl_cpardiso->n = A->rmap->N;
499:   mat_mkl_cpardiso->msglvl = 0;
500:   mat_mkl_cpardiso->nrhs = 1;
501:   mat_mkl_cpardiso->err = 0;
502:   mat_mkl_cpardiso->phase = -1;
503: #if defined(PETSC_USE_COMPLEX)
504:   mat_mkl_cpardiso->mtype = 13;
505: #else
506:   mat_mkl_cpardiso->mtype = 11;
507: #endif

509: #if defined(PETSC_USE_REAL_SINGLE)
510:   mat_mkl_cpardiso->iparm[27] = 1;
511: #else
512:   mat_mkl_cpardiso->iparm[27] = 0;
513: #endif

515:   mat_mkl_cpardiso->iparm[34] = 1;  /* C style */

517:   mat_mkl_cpardiso->iparm[ 0] =  1; /* Solver default parameters overriden with provided by iparm */
518:   mat_mkl_cpardiso->iparm[ 1] =  2; /* Use METIS for fill-in reordering */
519:   mat_mkl_cpardiso->iparm[ 5] =  0; /* Write solution into x */
520:   mat_mkl_cpardiso->iparm[ 7] =  2; /* Max number of iterative refinement steps */
521:   mat_mkl_cpardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
522:   mat_mkl_cpardiso->iparm[10] =  1; /* Use nonsymmetric permutation and scaling MPS */
523:   mat_mkl_cpardiso->iparm[12] =  1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
524:   mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
525:   mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
526:   mat_mkl_cpardiso->iparm[26] =  1; /* Check input data for correctness */

528:   mat_mkl_cpardiso->iparm[39] = 0;
529:   if (size > 1) {
530:     mat_mkl_cpardiso->iparm[39] = 2;
531:     mat_mkl_cpardiso->iparm[40] = A->rmap->rstart;
532:     mat_mkl_cpardiso->iparm[41] = A->rmap->rend-1;
533:   }
534:   mat_mkl_cpardiso->perm = 0;
535:   return(0);
536: }

538: /*
539:  * Symbolic decomposition. Mkl_Pardiso analysis phase.
540:  */
541: PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
542: {
543:   Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
544:   PetscErrorCode  ierr;

547:   mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;

549:   /* Set MKL_CPARDISO options from the options database */
550:   PetscSetMKL_CPARDISOFromOptions(F,A);
551:   (*mat_mkl_cpardiso->ConvertToTriples)(A,MAT_INITIAL_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);

553:   mat_mkl_cpardiso->n = A->rmap->N;

555:   /* analysis phase */
556:   /*----------------*/
557:   mat_mkl_cpardiso->phase = JOB_ANALYSIS;

559:   cluster_sparse_solver (
560:     mat_mkl_cpardiso->pt,
561:     &mat_mkl_cpardiso->maxfct,
562:     &mat_mkl_cpardiso->mnum,
563:     &mat_mkl_cpardiso->mtype,
564:     &mat_mkl_cpardiso->phase,
565:     &mat_mkl_cpardiso->n,
566:     mat_mkl_cpardiso->a,
567:     mat_mkl_cpardiso->ia,
568:     mat_mkl_cpardiso->ja,
569:     mat_mkl_cpardiso->perm,
570:     &mat_mkl_cpardiso->nrhs,
571:     mat_mkl_cpardiso->iparm,
572:     &mat_mkl_cpardiso->msglvl,
573:     NULL,
574:     NULL,
575:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
576:     (PetscInt*)&mat_mkl_cpardiso->err);

578:   if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\".Check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));

580:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
581:   F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO;
582:   F->ops->solve           = MatSolve_MKL_CPARDISO;
583:   F->ops->solvetranspose  = MatSolveTranspose_MKL_CPARDISO;
584:   F->ops->matsolve        = MatMatSolve_MKL_CPARDISO;
585:   return(0);
586: }

588: PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer)
589: {
590:   PetscErrorCode    ierr;
591:   PetscBool         iascii;
592:   PetscViewerFormat format;
593:   Mat_MKL_CPARDISO  *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
594:   PetscInt          i;

597:   /* check if matrix is mkl_cpardiso type */
598:   if (A->ops->solve != MatSolve_MKL_CPARDISO) return(0);

600:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
601:   if (iascii) {
602:     PetscViewerGetFormat(viewer,&format);
603:     if (format == PETSC_VIEWER_ASCII_INFO) {
604:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO run parameters:\n");
605:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO phase:             %d \n",mat_mkl_cpardiso->phase);
606:       for(i = 1; i <= 64; i++){
607:         PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO iparm[%d]:     %d \n",i, mat_mkl_cpardiso->iparm[i - 1]);
608:       }
609:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO maxfct:     %d \n", mat_mkl_cpardiso->maxfct);
610:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mnum:     %d \n", mat_mkl_cpardiso->mnum);
611:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mtype:     %d \n", mat_mkl_cpardiso->mtype);
612:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO n:     %d \n", mat_mkl_cpardiso->n);
613:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO nrhs:     %d \n", mat_mkl_cpardiso->nrhs);
614:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO msglvl:     %d \n", mat_mkl_cpardiso->msglvl);
615:     }
616:   }
617:   return(0);
618: }

620: PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info)
621: {
622:   Mat_MKL_CPARDISO *mat_mkl_cpardiso =(Mat_MKL_CPARDISO*)A->data;

625:   info->block_size        = 1.0;
626:   info->nz_allocated      = mat_mkl_cpardiso->nz + 0.0;
627:   info->nz_unneeded       = 0.0;
628:   info->assemblies        = 0.0;
629:   info->mallocs           = 0.0;
630:   info->memory            = 0.0;
631:   info->fill_ratio_given  = 0;
632:   info->fill_ratio_needed = 0;
633:   info->factor_mallocs    = 0;
634:   return(0);
635: }

637: PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F,PetscInt icntl,PetscInt ival)
638: {
639:   Mat_MKL_CPARDISO *mat_mkl_cpardiso =(Mat_MKL_CPARDISO*)F->data;

642:   if(icntl <= 64){
643:     mat_mkl_cpardiso->iparm[icntl - 1] = ival;
644:   } else {
645:     if(icntl == 65) mkl_set_num_threads((int)ival);
646:     else if(icntl == 66) mat_mkl_cpardiso->maxfct = ival;
647:     else if(icntl == 67) mat_mkl_cpardiso->mnum = ival;
648:     else if(icntl == 68) mat_mkl_cpardiso->msglvl = ival;
649:     else if(icntl == 69){
650:       mat_mkl_cpardiso->mtype = ival;
651: #if defined(PETSC_USE_REAL_SINGLE)
652:       mat_mkl_cpardiso->iparm[27] = 1;
653: #else
654:       mat_mkl_cpardiso->iparm[27] = 0;
655: #endif
656:       mat_mkl_cpardiso->iparm[34] = 1;
657:     }
658:   }
659:   return(0);
660: }

662: /*@
663:   MatMkl_CPardisoSetCntl - Set Mkl_Pardiso parameters

665:    Logically Collective on Mat

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

672:   Options Database:
673: .   -mat_mkl_cpardiso_<icntl> <ival>

675:    Level: Intermediate

677:    Notes:
678:     This routine cannot be used if you are solving the linear system with TS, SNES, or KSP, only if you directly call MatGetFactor() so use the options 
679:           database approach when working with TS, SNES, or KSP.

681:    References:
682: .      Mkl_Pardiso Users' Guide

684: .seealso: MatGetFactor()
685: @*/
686: PetscErrorCode MatMkl_CPardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
687: {

691:   PetscTryMethod(F,"MatMkl_CPardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
692:   return(0);
693: }

695: static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type)
696: {
698:   *type = MATSOLVERMKL_CPARDISO;
699:   return(0);
700: }

702: /* MatGetFactor for MPI AIJ matrices */
703: static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A,MatFactorType ftype,Mat *F)
704: {
705:   Mat              B;
706:   PetscErrorCode   ierr;
707:   Mat_MKL_CPARDISO *mat_mkl_cpardiso;
708:   PetscBool        isSeqAIJ;

711:   /* Create the factorization matrix */

713:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
714:   MatCreate(PetscObjectComm((PetscObject)A),&B);
715:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
716:   PetscStrallocpy("mkl_cpardiso",&((PetscObject)B)->type_name);
717:   MatSetUp(B);

719:   PetscNewLog(B,&mat_mkl_cpardiso);

721:   if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO;
722:   else          mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO;

724:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO;
725:   B->ops->destroy = MatDestroy_MKL_CPARDISO;

727:   B->ops->view    = MatView_MKL_CPARDISO;
728:   B->ops->getinfo = MatGetInfo_MKL_CPARDISO;

730:   B->factortype   = ftype;
731:   B->assembled    = PETSC_TRUE;           /* required by -ksp_view */

733:   B->data = mat_mkl_cpardiso;

735:   /* set solvertype */
736:   PetscFree(B->solvertype);
737:   PetscStrallocpy(MATSOLVERMKL_CPARDISO,&B->solvertype);

739:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_cpardiso);
740:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_CPardisoSetCntl_C",MatMkl_CPardisoSetCntl_MKL_CPARDISO);
741:   PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso);

743:   *F = B;
744:   return(0);
745: }

747: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void)
748: {
750: 
752:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
753:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
754:   return(0);
755: }