Actual source code: mkl_cpardiso.c

petsc-3.10.3 2018-12-18
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/aij/mpi/mpiaij.h>

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

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

 33: #define IPARM_SIZE 64
 34: #define INT_TYPE MKL_INT

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

 65: /*
 66:  *  Internal data structure.
 67:  *  For more information check mkl_cpardiso manual.
 68:  */

 70: typedef struct {

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

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

 81:   MPI_Comm     comm_mkl_cpardiso;

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

 86:   /* Matrix structure */
 87:   PetscScalar  *a;

 89:   INT_TYPE     *ia, *ja;

 91:   /* Number of non-zero elements */
 92:   INT_TYPE     nz;

 94:   /* Row permutaton vector*/
 95:   INT_TYPE     *perm;

 97:   /* Define is matrix preserve sparce structure. */
 98:   MatStructure matstruc;

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

102:   /* True if mkl_cpardiso function have been used. */
103:   PetscBool CleanUp;
104: } Mat_MKL_CPARDISO;

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

127:   *v=aa->a;
128:   if (reuse == MAT_INITIAL_MATRIX) {
129:     *r   = (INT_TYPE*)aa->i;
130:     *c   = (INT_TYPE*)aa->j;
131:     *nnz = aa->nz;
132:   }
133:   return(0);
134: }

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

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

153:   garray = mat->garray;

155:   if (reuse == MAT_INITIAL_MATRIX) {
156:     nz   = aa->nz + bb->nz;
157:     *nnz = nz;
158:     PetscMalloc((nz*(sizeof(PetscInt)+sizeof(PetscScalar)) + (m+1)*sizeof(PetscInt)), &row);
159:     col  = row + m + 1;
160:     val  = (PetscScalar*)(col + nz);
161:     *r = row; *c = col; *v = val;
162:     row[0] = 0;
163:   } else {
164:     row = *r; col = *c; val = *v;
165:   }

167:   nz = 0;
168:   for (i=0; i<m; i++) {
169:     row[i] = nz;
170:     countA     = ai[i+1] - ai[i];
171:     countB     = bi[i+1] - bi[i];
172:     ajj        = aj + ai[i]; /* ptr to the beginning of this row */
173:     bjj        = bj + bi[i];

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

189:     /* A part */
190:     for (j=0; j<countA; j++) {
191:       col[nz]   = rstart + ajj[j];
192:       val[nz++] = *av++;
193:     }

195:     /* B part, larger col index */
196:     for (j=jB; j<countB; j++) {
197:       col[nz]   = garray[bjj[j]];
198:       val[nz++] = *bv++;
199:     }
200:   }
201:   row[m] = nz;

203:   return(0);
204: }

206: /*
207:  * Free memory for Mat_MKL_CPARDISO structure and pointers to objects.
208:  */
209: PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A)
210: {
211:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
212:   PetscErrorCode   ierr;

215:   /* Terminate instance, deallocate memories */
216:   if (mat_mkl_cpardiso->CleanUp) {
217:     mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

219:     cluster_sparse_solver (
220:       mat_mkl_cpardiso->pt,
221:       &mat_mkl_cpardiso->maxfct,
222:       &mat_mkl_cpardiso->mnum,
223:       &mat_mkl_cpardiso->mtype,
224:       &mat_mkl_cpardiso->phase,
225:       &mat_mkl_cpardiso->n,
226:       NULL,
227:       NULL,
228:       NULL,
229:       mat_mkl_cpardiso->perm,
230:       &mat_mkl_cpardiso->nrhs,
231:       mat_mkl_cpardiso->iparm,
232:       &mat_mkl_cpardiso->msglvl,
233:       NULL,
234:       NULL,
235:       &mat_mkl_cpardiso->comm_mkl_cpardiso,
236:       (int*)&mat_mkl_cpardiso->err);
237:   }

239:   if (mat_mkl_cpardiso->ConvertToTriples == MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO) {
240:     PetscFree(mat_mkl_cpardiso->ia);
241:   }
242:   MPI_Comm_free(&(mat_mkl_cpardiso->comm_mkl_cpardiso));
243:   PetscFree(A->data);

245:   /* clear composed functions */
246:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
247:   PetscObjectComposeFunction((PetscObject)A,"MatMkl_CPardisoSetCntl_C",NULL);
248:   return(0);
249: }

251: /*
252:  * Computes Ax = b
253:  */
254: PetscErrorCode MatSolve_MKL_CPARDISO(Mat A,Vec b,Vec x)
255: {
256:   Mat_MKL_CPARDISO   *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
257:   PetscErrorCode    ierr;
258:   PetscScalar       *xarray;
259:   const PetscScalar *barray;

262:   mat_mkl_cpardiso->nrhs = 1;
263:   VecGetArray(x,&xarray);
264:   VecGetArrayRead(b,&barray);

266:   /* solve phase */
267:   /*-------------*/
268:   mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
269:   cluster_sparse_solver (
270:     mat_mkl_cpardiso->pt,
271:     &mat_mkl_cpardiso->maxfct,
272:     &mat_mkl_cpardiso->mnum,
273:     &mat_mkl_cpardiso->mtype,
274:     &mat_mkl_cpardiso->phase,
275:     &mat_mkl_cpardiso->n,
276:     mat_mkl_cpardiso->a,
277:     mat_mkl_cpardiso->ia,
278:     mat_mkl_cpardiso->ja,
279:     mat_mkl_cpardiso->perm,
280:     &mat_mkl_cpardiso->nrhs,
281:     mat_mkl_cpardiso->iparm,
282:     &mat_mkl_cpardiso->msglvl,
283:     (void*)barray,
284:     (void*)xarray,
285:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
286:     (int*)&mat_mkl_cpardiso->err);

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

290:   VecRestoreArray(x,&xarray);
291:   VecRestoreArrayRead(b,&barray);
292:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
293:   return(0);
294: }

296: PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A,Vec b,Vec x)
297: {
298:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
299:   PetscErrorCode   ierr;

302: #if defined(PETSC_USE_COMPLEX)
303:   mat_mkl_cpardiso->iparm[12 - 1] = 1;
304: #else
305:   mat_mkl_cpardiso->iparm[12 - 1] = 2;
306: #endif
307:   MatSolve_MKL_CPARDISO(A,b,x);
308:   mat_mkl_cpardiso->iparm[12 - 1] = 0;
309:   return(0);
310: }

312: PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A,Mat B,Mat X)
313: {
314:   Mat_MKL_CPARDISO  *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
315:   PetscErrorCode    ierr;
316:   PetscScalar       *xarray;
317:   const PetscScalar *barray;

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

322:   if(mat_mkl_cpardiso->nrhs > 0){
323:     MatDenseGetArrayRead(B,&barray);
324:     MatDenseGetArray(X,&xarray);

326:     /* solve phase */
327:     /*-------------*/
328:     mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
329:     cluster_sparse_solver (
330:       mat_mkl_cpardiso->pt,
331:       &mat_mkl_cpardiso->maxfct,
332:       &mat_mkl_cpardiso->mnum,
333:       &mat_mkl_cpardiso->mtype,
334:       &mat_mkl_cpardiso->phase,
335:       &mat_mkl_cpardiso->n,
336:       mat_mkl_cpardiso->a,
337:       mat_mkl_cpardiso->ia,
338:       mat_mkl_cpardiso->ja,
339:       mat_mkl_cpardiso->perm,
340:       &mat_mkl_cpardiso->nrhs,
341:       mat_mkl_cpardiso->iparm,
342:       &mat_mkl_cpardiso->msglvl,
343:       (void*)barray,
344:       (void*)xarray,
345:       &mat_mkl_cpardiso->comm_mkl_cpardiso,
346:       (int*)&mat_mkl_cpardiso->err);
347:     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));
348:     MatDenseRestoreArrayRead(B,&barray);
349:     MatDenseRestoreArray(X,&xarray);

351:   }
352:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
353:   return(0);

355: }

357: /*
358:  * LU Decomposition
359:  */
360: PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F,Mat A,const MatFactorInfo *info)
361: {
362:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(F)->data;
363:   PetscErrorCode   ierr;

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

369:   mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION;
370:   cluster_sparse_solver (
371:     mat_mkl_cpardiso->pt,
372:     &mat_mkl_cpardiso->maxfct,
373:     &mat_mkl_cpardiso->mnum,
374:     &mat_mkl_cpardiso->mtype,
375:     &mat_mkl_cpardiso->phase,
376:     &mat_mkl_cpardiso->n,
377:     mat_mkl_cpardiso->a,
378:     mat_mkl_cpardiso->ia,
379:     mat_mkl_cpardiso->ja,
380:     mat_mkl_cpardiso->perm,
381:     &mat_mkl_cpardiso->nrhs,
382:     mat_mkl_cpardiso->iparm,
383:     &mat_mkl_cpardiso->msglvl,
384:     NULL,
385:     NULL,
386:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
387:     &mat_mkl_cpardiso->err);
388:   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));

390:   mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
391:   mat_mkl_cpardiso->CleanUp  = PETSC_TRUE;
392:   return(0);
393: }

395: /* Sets mkl_cpardiso options from the options database */
396: PetscErrorCode PetscSetMKL_CPARDISOFromOptions(Mat F, Mat A)
397: {
398:   Mat_MKL_CPARDISO    *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
399:   PetscErrorCode      ierr;
400:   PetscInt            icntl;
401:   PetscBool           flg;
402:   int                 threads;

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

409:   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);
410:   if (flg) mat_mkl_cpardiso->maxfct = icntl;

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

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

418:   PetscOptionsInt("-mat_mkl_cpardiso_69","Defines the matrix type","None",mat_mkl_cpardiso->mtype,&icntl,&flg);
419:   if(flg){
420:     mat_mkl_cpardiso->mtype = icntl;
421: #if defined(PETSC_USE_REAL_SINGLE)
422:     mat_mkl_cpardiso->iparm[27] = 1;
423: #else
424:     mat_mkl_cpardiso->iparm[27] = 0;
425: #endif
426:     mat_mkl_cpardiso->iparm[34] = 1;
427:   }
428:   PetscOptionsInt("-mat_mkl_cpardiso_1","Use default values","None",mat_mkl_cpardiso->iparm[0],&icntl,&flg);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

488:   PetscOptionsEnd();
489:   return(0);
490: }

492: PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso)
493: {
494:   PetscErrorCode  ierr;
495:   PetscMPIInt     size;


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

502:   mat_mkl_cpardiso->CleanUp = PETSC_FALSE;
503:   mat_mkl_cpardiso->maxfct = 1;
504:   mat_mkl_cpardiso->mnum = 1;
505:   mat_mkl_cpardiso->n = A->rmap->N;
506:   mat_mkl_cpardiso->msglvl = 0;
507:   mat_mkl_cpardiso->nrhs = 1;
508:   mat_mkl_cpardiso->err = 0;
509:   mat_mkl_cpardiso->phase = -1;
510: #if defined(PETSC_USE_COMPLEX)
511:   mat_mkl_cpardiso->mtype = 13;
512: #else
513:   mat_mkl_cpardiso->mtype = 11;
514: #endif

516: #if defined(PETSC_USE_REAL_SINGLE)
517:   mat_mkl_cpardiso->iparm[27] = 1;
518: #else
519:   mat_mkl_cpardiso->iparm[27] = 0;
520: #endif

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

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

535:   mat_mkl_cpardiso->iparm[39] = 0;
536:   if (size > 1) {
537:     mat_mkl_cpardiso->iparm[39] = 2;
538:     mat_mkl_cpardiso->iparm[40] = A->rmap->rstart;
539:     mat_mkl_cpardiso->iparm[41] = A->rmap->rend-1;
540:   }
541:   mat_mkl_cpardiso->perm = 0;
542:   return(0);
543: }

545: /*
546:  * Symbolic decomposition. Mkl_Pardiso analysis phase.
547:  */
548: PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
549: {
550:   Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
551:   PetscErrorCode  ierr;

554:   mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;

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

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

562:   /* analysis phase */
563:   /*----------------*/
564:   mat_mkl_cpardiso->phase = JOB_ANALYSIS;

566:   cluster_sparse_solver (
567:     mat_mkl_cpardiso->pt,
568:     &mat_mkl_cpardiso->maxfct,
569:     &mat_mkl_cpardiso->mnum,
570:     &mat_mkl_cpardiso->mtype,
571:     &mat_mkl_cpardiso->phase,
572:     &mat_mkl_cpardiso->n,
573:     mat_mkl_cpardiso->a,
574:     mat_mkl_cpardiso->ia,
575:     mat_mkl_cpardiso->ja,
576:     mat_mkl_cpardiso->perm,
577:     &mat_mkl_cpardiso->nrhs,
578:     mat_mkl_cpardiso->iparm,
579:     &mat_mkl_cpardiso->msglvl,
580:     NULL,
581:     NULL,
582:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
583:     (int*)&mat_mkl_cpardiso->err);

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

587:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
588:   F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO;
589:   F->ops->solve           = MatSolve_MKL_CPARDISO;
590:   F->ops->solvetranspose  = MatSolveTranspose_MKL_CPARDISO;
591:   F->ops->matsolve        = MatMatSolve_MKL_CPARDISO;
592:   return(0);
593: }

595: PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer)
596: {
597:   PetscErrorCode    ierr;
598:   PetscBool         iascii;
599:   PetscViewerFormat format;
600:   Mat_MKL_CPARDISO  *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
601:   PetscInt          i;

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

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

627: PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info)
628: {
629:   Mat_MKL_CPARDISO *mat_mkl_cpardiso =(Mat_MKL_CPARDISO*)A->data;

632:   info->block_size        = 1.0;
633:   info->nz_allocated      = mat_mkl_cpardiso->nz + 0.0;
634:   info->nz_unneeded       = 0.0;
635:   info->assemblies        = 0.0;
636:   info->mallocs           = 0.0;
637:   info->memory            = 0.0;
638:   info->fill_ratio_given  = 0;
639:   info->fill_ratio_needed = 0;
640:   info->factor_mallocs    = 0;
641:   return(0);
642: }

644: PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F,PetscInt icntl,PetscInt ival)
645: {
646:   Mat_MKL_CPARDISO *mat_mkl_cpardiso =(Mat_MKL_CPARDISO*)F->data;

649:   if(icntl <= 64){
650:     mat_mkl_cpardiso->iparm[icntl - 1] = ival;
651:   } else {
652:     if(icntl == 65) mkl_set_num_threads((int)ival);
653:     else if(icntl == 66) mat_mkl_cpardiso->maxfct = ival;
654:     else if(icntl == 67) mat_mkl_cpardiso->mnum = ival;
655:     else if(icntl == 68) mat_mkl_cpardiso->msglvl = ival;
656:     else if(icntl == 69){
657:       mat_mkl_cpardiso->mtype = ival;
658: #if defined(PETSC_USE_REAL_SINGLE)
659:       mat_mkl_cpardiso->iparm[27] = 1;
660: #else
661:       mat_mkl_cpardiso->iparm[27] = 0;
662: #endif
663:       mat_mkl_cpardiso->iparm[34] = 1;
664:     }
665:   }
666:   return(0);
667: }

669: /*@
670:   MatMkl_CPardisoSetCntl - Set Mkl_Pardiso parameters

672:    Logically Collective on Mat

674:    Input Parameters:
675: +  F - the factored matrix obtained by calling MatGetFactor()
676: .  icntl - index of Mkl_Pardiso parameter
677: -  ival - value of Mkl_Pardiso parameter

679:   Options Database:
680: .   -mat_mkl_cpardiso_<icntl> <ival>

682:    Level: Intermediate

684:    Notes:
685:     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 
686:           database approach when working with TS, SNES, or KSP.

688:    References:
689: .      Mkl_Pardiso Users' Guide

691: .seealso: MatGetFactor()
692: @*/
693: PetscErrorCode MatMkl_CPardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
694: {

698:   PetscTryMethod(F,"MatMkl_CPardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
699:   return(0);
700: }

702: static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type)
703: {
705:   *type = MATSOLVERMKL_CPARDISO;
706:   return(0);
707: }

709: /* MatGetFactor for MPI AIJ matrices */
710: static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A,MatFactorType ftype,Mat *F)
711: {
712:   Mat              B;
713:   PetscErrorCode   ierr;
714:   Mat_MKL_CPARDISO *mat_mkl_cpardiso;
715:   PetscBool        isSeqAIJ;

718:   /* Create the factorization matrix */

720:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
721:   MatCreate(PetscObjectComm((PetscObject)A),&B);
722:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
723:   PetscStrallocpy("mkl_cpardiso",&((PetscObject)B)->type_name);
724:   MatSetUp(B);

726:   PetscNewLog(B,&mat_mkl_cpardiso);

728:   if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO;
729:   else          mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO;

731:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO;
732:   B->ops->destroy = MatDestroy_MKL_CPARDISO;

734:   B->ops->view    = MatView_MKL_CPARDISO;
735:   B->ops->getinfo = MatGetInfo_MKL_CPARDISO;

737:   B->factortype   = ftype;
738:   B->assembled    = PETSC_TRUE;           /* required by -ksp_view */

740:   B->data = mat_mkl_cpardiso;

742:   /* set solvertype */
743:   PetscFree(B->solvertype);
744:   PetscStrallocpy(MATSOLVERMKL_CPARDISO,&B->solvertype);

746:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_cpardiso);
747:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_CPardisoSetCntl_C",MatMkl_CPardisoSetCntl_MKL_CPARDISO);
748:   PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso);

750:   *F = B;
751:   return(0);
752: }

754: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void)
755: {
757: 
759:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
760:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
761:   return(0);
762: }