Actual source code: mkl_pardiso.c
petsc-3.5.4 2015-05-23
1: #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
2: #include <../src/mat/impls/dense/seq/dense.h>
4: #include <stdio.h>
5: #include <stdlib.h>
6: #include <math.h>
7: #include <mkl.h>
9: /*
10: * Possible mkl_pardiso phases that controls the execution of the solver.
11: * For more information check mkl_pardiso manual.
12: */
13: #define JOB_ANALYSIS 11
14: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
15: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
16: #define JOB_NUMERICAL_FACTORIZATION 22
17: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
18: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
19: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
20: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
21: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
22: #define JOB_RELEASE_OF_LU_MEMORY 0
23: #define JOB_RELEASE_OF_ALL_MEMORY -1
25: #define IPARM_SIZE 64
26: #if defined(PETSC_USE_64BIT_INDICES)
27: #define INT_TYPE long long int
28: #define MKL_PARDISO pardiso_64
29: #define MKL_PARDISO_INIT pardiso_64init
30: #else
31: #define INT_TYPE int
32: #define MKL_PARDISO pardiso
33: #define MKL_PARDISO_INIT pardisoinit
34: #endif
37: /*
38: * Internal data structure.
39: * For more information check mkl_pardiso manual.
40: */
41: typedef struct {
43: /*Configuration vector*/
44: INT_TYPE iparm[IPARM_SIZE];
46: /*
47: * Internal mkl_pardiso memory location.
48: * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
49: */
50: void *pt[IPARM_SIZE];
52: /*Basic mkl_pardiso info*/
53: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
55: /*Matrix structure*/
56: void *a;
57: INT_TYPE *ia, *ja;
59: /*Number of non-zero elements*/
60: INT_TYPE nz;
62: /*Row permutaton vector*/
63: INT_TYPE *perm;
65: /*Deffine is matrix preserve sparce structure.*/
66: MatStructure matstruc;
68: /*True if mkl_pardiso function have been used.*/
69: PetscBool CleanUp;
70: } Mat_MKL_PARDISO;
73: void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm []){
74: int iparm_copy[IPARM_SIZE], mtype_copy, i;
75: mtype_copy = *mtype;
76: pardisoinit(pt, &mtype_copy, iparm_copy);
77: for(i = 0; i < IPARM_SIZE; i++)
78: iparm[i] = iparm_copy[i];
79: }
82: /*
83: * Copy the elements of matrix A.
84: * Input:
85: * - Mat A: MATSEQAIJ matrix
86: * - int shift: matrix index.
87: * - 0 for c representation
88: * - 1 for fortran representation
89: * - MatReuse reuse:
90: * - MAT_INITIAL_MATRIX: Create a new aij representation
91: * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
92: * Output:
93: * - int *nnz: Number of nonzero-elements.
94: * - int **r pointer to i index
95: * - int **c pointer to j elements
96: * - MATRIXTYPE **v: Non-zero elements
97: */
100: PetscErrorCode MatCopy_MKL_PARDISO(Mat A, MatReuse reuse, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, void **v){
102: Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data;
105: *v=aa->a;
106: if (reuse == MAT_INITIAL_MATRIX) {
107: *r = (INT_TYPE*)aa->i;
108: *c = (INT_TYPE*)aa->j;
109: *nnz = aa->nz;
110: }
111: return(0);
112: }
115: /*
116: * Free memory for Mat_MKL_PARDISO structure and pointers to objects.
117: */
120: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A){
121: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
125: /* Terminate instance, deallocate memories */
126: if (mat_mkl_pardiso->CleanUp) {
127: mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
130: MKL_PARDISO (mat_mkl_pardiso->pt,
131: &mat_mkl_pardiso->maxfct,
132: &mat_mkl_pardiso->mnum,
133: &mat_mkl_pardiso->mtype,
134: &mat_mkl_pardiso->phase,
135: &mat_mkl_pardiso->n,
136: NULL,
137: NULL,
138: NULL,
139: mat_mkl_pardiso->perm,
140: &mat_mkl_pardiso->nrhs,
141: mat_mkl_pardiso->iparm,
142: &mat_mkl_pardiso->msglvl,
143: NULL,
144: NULL,
145: &mat_mkl_pardiso->err);
146: }
147: PetscFree(mat_mkl_pardiso->perm);
148: PetscFree(A->spptr);
150: /* clear composed functions */
151: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
152: PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);
153: return(0);
154: }
157: /*
158: * Computes Ax = b
159: */
162: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x){
163: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
164: PetscErrorCode ierr;
165: PetscScalar *barray, *xarray;
170: mat_mkl_pardiso->nrhs = 1;
171: VecGetArray(x,&xarray);
172: VecGetArray(b,&barray);
174: /* solve phase */
175: /*-------------*/
176: mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
177: MKL_PARDISO (mat_mkl_pardiso->pt,
178: &mat_mkl_pardiso->maxfct,
179: &mat_mkl_pardiso->mnum,
180: &mat_mkl_pardiso->mtype,
181: &mat_mkl_pardiso->phase,
182: &mat_mkl_pardiso->n,
183: mat_mkl_pardiso->a,
184: mat_mkl_pardiso->ia,
185: mat_mkl_pardiso->ja,
186: mat_mkl_pardiso->perm,
187: &mat_mkl_pardiso->nrhs,
188: mat_mkl_pardiso->iparm,
189: &mat_mkl_pardiso->msglvl,
190: (void*)barray,
191: (void*)xarray,
192: &mat_mkl_pardiso->err);
195: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err);
197: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
198: return(0);
199: }
204: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x){
205: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
209: #if defined(PETSC_USE_COMPLEX)
210: mat_mkl_pardiso->iparm[12 - 1] = 1;
211: #else
212: mat_mkl_pardiso->iparm[12 - 1] = 2;
213: #endif
214: MatSolve_MKL_PARDISO(A,b,x);
215: mat_mkl_pardiso->iparm[12 - 1] = 0;
216: return(0);
217: }
222: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X){
223: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
224: PetscErrorCode ierr;
225: PetscScalar *barray, *xarray;
226: PetscBool flg;
230: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
231: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
232: PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
233: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");
235: MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);
237: if(mat_mkl_pardiso->nrhs > 0){
238: MatDenseGetArray(B,&barray);
239: MatDenseGetArray(X,&xarray);
241: /* solve phase */
242: /*-------------*/
243: mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
244: MKL_PARDISO (mat_mkl_pardiso->pt,
245: &mat_mkl_pardiso->maxfct,
246: &mat_mkl_pardiso->mnum,
247: &mat_mkl_pardiso->mtype,
248: &mat_mkl_pardiso->phase,
249: &mat_mkl_pardiso->n,
250: mat_mkl_pardiso->a,
251: mat_mkl_pardiso->ia,
252: mat_mkl_pardiso->ja,
253: mat_mkl_pardiso->perm,
254: &mat_mkl_pardiso->nrhs,
255: mat_mkl_pardiso->iparm,
256: &mat_mkl_pardiso->msglvl,
257: (void*)barray,
258: (void*)xarray,
259: &mat_mkl_pardiso->err);
260: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err);
261: }
262: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
263: return(0);
265: }
267: /*
268: * LU Decomposition
269: */
272: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info){
273: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->spptr;
276: /* numerical factorization phase */
277: /*-------------------------------*/
281: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
282: MatCopy_MKL_PARDISO(A, MAT_REUSE_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);
284: /* numerical factorization phase */
285: /*-------------------------------*/
286: mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
287: MKL_PARDISO (mat_mkl_pardiso->pt,
288: &mat_mkl_pardiso->maxfct,
289: &mat_mkl_pardiso->mnum,
290: &mat_mkl_pardiso->mtype,
291: &mat_mkl_pardiso->phase,
292: &mat_mkl_pardiso->n,
293: mat_mkl_pardiso->a,
294: mat_mkl_pardiso->ia,
295: mat_mkl_pardiso->ja,
296: mat_mkl_pardiso->perm,
297: &mat_mkl_pardiso->nrhs,
298: mat_mkl_pardiso->iparm,
299: &mat_mkl_pardiso->msglvl,
300: NULL,
301: NULL,
302: &mat_mkl_pardiso->err);
303: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err);
305: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
306: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
307: return(0);
308: }
310: /* Sets mkl_pardiso options from the options database */
313: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A){
314: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
315: PetscErrorCode ierr;
316: PetscInt icntl;
317: PetscBool flg;
318: int pt[IPARM_SIZE], threads;
321: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");
322: PetscOptionsInt("-mat_mkl_pardiso_65",
323: "Number of thrads to use",
324: "None",
325: threads,
326: &threads,
327: &flg);
328: if (flg) mkl_set_num_threads(threads);
330: PetscOptionsInt("-mat_mkl_pardiso_66",
331: "Maximum number of factors with identical sparsity structure that must be kept in memory at the same time",
332: "None",
333: mat_mkl_pardiso->maxfct,
334: &icntl,
335: &flg);
336: if (flg) mat_mkl_pardiso->maxfct = icntl;
338: PetscOptionsInt("-mat_mkl_pardiso_67",
339: "Indicates the actual matrix for the solution phase",
340: "None",
341: mat_mkl_pardiso->mnum,
342: &icntl,
343: &flg);
344: if (flg) mat_mkl_pardiso->mnum = icntl;
345:
346: PetscOptionsInt("-mat_mkl_pardiso_68",
347: "Message level information",
348: "None",
349: mat_mkl_pardiso->msglvl,
350: &icntl,
351: &flg);
352: if (flg) mat_mkl_pardiso->msglvl = icntl;
354: PetscOptionsInt("-mat_mkl_pardiso_69",
355: "Defines the matrix type",
356: "None",
357: mat_mkl_pardiso->mtype,
358: &icntl,
359: &flg);
360: if(flg){
361: mat_mkl_pardiso->mtype = icntl;
362: MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
363: #if defined(PETSC_USE_REAL_SINGLE)
364: mat_mkl_pardiso->iparm[27] = 1;
365: #else
366: mat_mkl_pardiso->iparm[27] = 0;
367: #endif
368: mat_mkl_pardiso->iparm[34] = 1;
369: }
370: PetscOptionsInt("-mat_mkl_pardiso_1",
371: "Use default values",
372: "None",
373: mat_mkl_pardiso->iparm[0],
374: &icntl,
375: &flg);
377: if(flg && icntl != 0){
378: PetscOptionsInt("-mat_mkl_pardiso_2",
379: "Fill-in reducing ordering for the input matrix",
380: "None",
381: mat_mkl_pardiso->iparm[1],
382: &icntl,
383: &flg);
384: if (flg) mat_mkl_pardiso->iparm[1] = icntl;
386: PetscOptionsInt("-mat_mkl_pardiso_4",
387: "Preconditioned CGS/CG",
388: "None",
389: mat_mkl_pardiso->iparm[3],
390: &icntl,
391: &flg);
392: if (flg) mat_mkl_pardiso->iparm[3] = icntl;
394: PetscOptionsInt("-mat_mkl_pardiso_5",
395: "User permutation",
396: "None",
397: mat_mkl_pardiso->iparm[4],
398: &icntl,
399: &flg);
400: if (flg) mat_mkl_pardiso->iparm[4] = icntl;
402: PetscOptionsInt("-mat_mkl_pardiso_6",
403: "Write solution on x",
404: "None",
405: mat_mkl_pardiso->iparm[5],
406: &icntl,
407: &flg);
408: if (flg) mat_mkl_pardiso->iparm[5] = icntl;
410: PetscOptionsInt("-mat_mkl_pardiso_8",
411: "Iterative refinement step",
412: "None",
413: mat_mkl_pardiso->iparm[7],
414: &icntl,
415: &flg);
416: if (flg) mat_mkl_pardiso->iparm[7] = icntl;
418: PetscOptionsInt("-mat_mkl_pardiso_10",
419: "Pivoting perturbation",
420: "None",
421: mat_mkl_pardiso->iparm[9],
422: &icntl,
423: &flg);
424: if (flg) mat_mkl_pardiso->iparm[9] = icntl;
426: PetscOptionsInt("-mat_mkl_pardiso_11",
427: "Scaling vectors",
428: "None",
429: mat_mkl_pardiso->iparm[10],
430: &icntl,
431: &flg);
432: if (flg) mat_mkl_pardiso->iparm[10] = icntl;
434: PetscOptionsInt("-mat_mkl_pardiso_12",
435: "Solve with transposed or conjugate transposed matrix A",
436: "None",
437: mat_mkl_pardiso->iparm[11],
438: &icntl,
439: &flg);
440: if (flg) mat_mkl_pardiso->iparm[11] = icntl;
442: PetscOptionsInt("-mat_mkl_pardiso_13",
443: "Improved accuracy using (non-) symmetric weighted matching",
444: "None",
445: mat_mkl_pardiso->iparm[12],
446: &icntl,
447: &flg);
448: if (flg) mat_mkl_pardiso->iparm[12] = icntl;
450: PetscOptionsInt("-mat_mkl_pardiso_18",
451: "Numbers of non-zero elements",
452: "None",
453: mat_mkl_pardiso->iparm[17],
454: &icntl,
455: &flg);
456: if (flg) mat_mkl_pardiso->iparm[17] = icntl;
458: PetscOptionsInt("-mat_mkl_pardiso_19",
459: "Report number of floating point operations",
460: "None",
461: mat_mkl_pardiso->iparm[18],
462: &icntl,
463: &flg);
464: if (flg) mat_mkl_pardiso->iparm[18] = icntl;
466: PetscOptionsInt("-mat_mkl_pardiso_21",
467: "Pivoting for symmetric indefinite matrices",
468: "None",
469: mat_mkl_pardiso->iparm[20],
470: &icntl,
471: &flg);
472: if (flg) mat_mkl_pardiso->iparm[20] = icntl;
474: PetscOptionsInt("-mat_mkl_pardiso_24",
475: "Parallel factorization control",
476: "None",
477: mat_mkl_pardiso->iparm[23],
478: &icntl,
479: &flg);
480: if (flg) mat_mkl_pardiso->iparm[23] = icntl;
482: PetscOptionsInt("-mat_mkl_pardiso_25",
483: "Parallel forward/backward solve control",
484: "None",
485: mat_mkl_pardiso->iparm[24],
486: &icntl,
487: &flg);
488: if (flg) mat_mkl_pardiso->iparm[24] = icntl;
490: PetscOptionsInt("-mat_mkl_pardiso_27",
491: "Matrix checker",
492: "None",
493: mat_mkl_pardiso->iparm[26],
494: &icntl,
495: &flg);
496: if (flg) mat_mkl_pardiso->iparm[26] = icntl;
498: PetscOptionsInt("-mat_mkl_pardiso_31",
499: "Partial solve and computing selected components of the solution vectors",
500: "None",
501: mat_mkl_pardiso->iparm[30],
502: &icntl,
503: &flg);
504: if (flg) mat_mkl_pardiso->iparm[30] = icntl;
506: PetscOptionsInt("-mat_mkl_pardiso_34",
507: "Optimal number of threads for conditional numerical reproducibility (CNR) mode",
508: "None",
509: mat_mkl_pardiso->iparm[33],
510: &icntl,
511: &flg);
512: if (flg) mat_mkl_pardiso->iparm[33] = icntl;
514: PetscOptionsInt("-mat_mkl_pardiso_60",
515: "Intel MKL_PARDISO mode",
516: "None",
517: mat_mkl_pardiso->iparm[59],
518: &icntl,
519: &flg);
520: if (flg) mat_mkl_pardiso->iparm[59] = icntl;
521: }
523: PetscOptionsEnd();
524: return(0);
525: }
530: PetscErrorCode PetscInitializeMKL_PARDISO(Mat A, Mat_MKL_PARDISO *mat_mkl_pardiso){
532: PetscInt i;
535: mat_mkl_pardiso->CleanUp = PETSC_FALSE;
536: mat_mkl_pardiso->maxfct = 1;
537: mat_mkl_pardiso->mnum = 1;
538: mat_mkl_pardiso->n = A->rmap->N;
539: mat_mkl_pardiso->msglvl = 0;
540: mat_mkl_pardiso->nrhs = 1;
541: mat_mkl_pardiso->err = 0;
542: mat_mkl_pardiso->phase = -1;
543: #if defined(PETSC_USE_COMPLEX)
544: mat_mkl_pardiso->mtype = 13;
545: #else
546: mat_mkl_pardiso->mtype = 11;
547: #endif
549: MKL_PARDISO_INIT(mat_mkl_pardiso->pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
551: #if defined(PETSC_USE_REAL_SINGLE)
552: mat_mkl_pardiso->iparm[27] = 1;
553: #else
554: mat_mkl_pardiso->iparm[27] = 0;
555: #endif
557: mat_mkl_pardiso->iparm[34] = 1;
559: PetscMalloc(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
560: for(i = 0; i < A->rmap->N; i++)
561: mat_mkl_pardiso->perm[i] = 0;
562: return(0);
563: }
566: /*
567: * Symbolic decomposition. Mkl_Pardiso analysis phase.
568: */
571: PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info){
573: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
578: mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
580: /* Set MKL_PARDISO options from the options database */
581: PetscSetMKL_PARDISOFromOptions(F,A);
583: MatCopy_MKL_PARDISO(A, MAT_INITIAL_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);
584: mat_mkl_pardiso->n = A->rmap->N;
586: /* analysis phase */
587: /*----------------*/
589: mat_mkl_pardiso->phase = JOB_ANALYSIS;
591: MKL_PARDISO (mat_mkl_pardiso->pt,
592: &mat_mkl_pardiso->maxfct,
593: &mat_mkl_pardiso->mnum,
594: &mat_mkl_pardiso->mtype,
595: &mat_mkl_pardiso->phase,
596: &mat_mkl_pardiso->n,
597: mat_mkl_pardiso->a,
598: mat_mkl_pardiso->ia,
599: mat_mkl_pardiso->ja,
600: mat_mkl_pardiso->perm,
601: &mat_mkl_pardiso->nrhs,
602: mat_mkl_pardiso->iparm,
603: &mat_mkl_pardiso->msglvl,
604: NULL,
605: NULL,
606: &mat_mkl_pardiso->err);
608: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d\n. Please check manual",mat_mkl_pardiso->err);
610: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
611: F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
612: F->ops->solve = MatSolve_MKL_PARDISO;
613: F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO;
614: F->ops->matsolve = MatMatSolve_MKL_PARDISO;
615: return(0);
616: }
621: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer){
622: PetscErrorCode ierr;
623: PetscBool iascii;
624: PetscViewerFormat format;
625: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
626: PetscInt i;
629: /* check if matrix is mkl_pardiso type */
630: if (A->ops->solve != MatSolve_MKL_PARDISO) return(0);
632: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
633: if (iascii) {
634: PetscViewerGetFormat(viewer,&format);
635: if (format == PETSC_VIEWER_ASCII_INFO) {
636: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
637: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);
638: for(i = 1; i <= 64; i++){
639: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
640: }
641: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);
642: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);
643: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);
644: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);
645: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);
646: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);
647: }
648: }
649: return(0);
650: }
655: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info){
656: Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->spptr;
659: info->block_size = 1.0;
660: info->nz_allocated = mat_mkl_pardiso->nz + 0.0;
661: info->nz_unneeded = 0.0;
662: info->assemblies = 0.0;
663: info->mallocs = 0.0;
664: info->memory = 0.0;
665: info->fill_ratio_given = 0;
666: info->fill_ratio_needed = 0;
667: info->factor_mallocs = 0;
668: return(0);
669: }
673: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival){
674: Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->spptr;
676: if(icntl <= 64){
677: mat_mkl_pardiso->iparm[icntl - 1] = ival;
678: } else {
679: if(icntl == 65)
680: mkl_set_num_threads((int)ival);
681: else if(icntl == 66)
682: mat_mkl_pardiso->maxfct = ival;
683: else if(icntl == 67)
684: mat_mkl_pardiso->mnum = ival;
685: else if(icntl == 68)
686: mat_mkl_pardiso->msglvl = ival;
687: else if(icntl == 69){
688: int pt[IPARM_SIZE];
689: mat_mkl_pardiso->mtype = ival;
690: MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
691: #if defined(PETSC_USE_REAL_SINGLE)
692: mat_mkl_pardiso->iparm[27] = 1;
693: #else
694: mat_mkl_pardiso->iparm[27] = 0;
695: #endif
696: mat_mkl_pardiso->iparm[34] = 1;
697: }
698: }
699: return(0);
700: }
704: /*@
705: MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters
707: Logically Collective on Mat
709: Input Parameters:
710: + F - the factored matrix obtained by calling MatGetFactor()
711: . icntl - index of Mkl_Pardiso parameter
712: - ival - value of Mkl_Pardiso parameter
714: Options Database:
715: . -mat_mkl_pardiso_<icntl> <ival>
717: Level: beginner
719: References: Mkl_Pardiso Users' Guide
721: .seealso: MatGetFactor()
722: @*/
723: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
724: {
728: PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
729: return(0);
730: }
733: /*MC
734: MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for
735: sequential matrices via the external package MKL_PARDISO.
737: Works with MATSEQAIJ matrices
739: Options Database Keys:
740: + -mat_mkl_pardiso_65 - Number of thrads to use
741: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
742: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
743: . -mat_mkl_pardiso_68 - Message level information
744: . -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
745: . -mat_mkl_pardiso_1 - Use default values
746: . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix
747: . -mat_mkl_pardiso_4 - Preconditioned CGS/CG
748: . -mat_mkl_pardiso_5 - User permutation
749: . -mat_mkl_pardiso_6 - Write solution on x
750: . -mat_mkl_pardiso_8 - Iterative refinement step
751: . -mat_mkl_pardiso_10 - Pivoting perturbation
752: . -mat_mkl_pardiso_11 - Scaling vectors
753: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
754: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
755: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
756: . -mat_mkl_pardiso_19 - Report number of floating point operations
757: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
758: . -mat_mkl_pardiso_24 - Parallel factorization control
759: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
760: . -mat_mkl_pardiso_27 - Matrix checker
761: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
762: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
763: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode
765: Level: beginner
767: For more information please check mkl_pardiso manual
769: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
771: M*/
776: static PetscErrorCode MatFactorGetSolverPackage_mkl_pardiso(Mat A, const MatSolverPackage *type){
778: *type = MATSOLVERMKL_PARDISO;
779: return(0);
780: }
783: /* MatGetFactor for Seq AIJ matrices */
786: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F){
787: Mat B;
789: Mat_MKL_PARDISO *mat_mkl_pardiso;
790: PetscBool isSeqAIJ;
793: /* Create the factorization matrix */
796: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
797: MatCreate(PetscObjectComm((PetscObject)A),&B);
798: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
799: MatSetType(B,((PetscObject)A)->type_name);
800: if (isSeqAIJ) {
801: MatSeqAIJSetPreallocation(B,0,NULL);
802: } else {
803: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Is not allowed other types of matrices apart from MATSEQAIJ.");
804: }
806: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
807: B->ops->destroy = MatDestroy_MKL_PARDISO;
808: B->ops->view = MatView_MKL_PARDISO;
809: B->factortype = ftype;
810: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
811: B->assembled = PETSC_TRUE; /* required by -ksp_view */
813: PetscNewLog(B,&mat_mkl_pardiso);
814: B->spptr = mat_mkl_pardiso;
816: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);
817: PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
818: PetscInitializeMKL_PARDISO(A, mat_mkl_pardiso);
820: *F = B;
821: return(0);
822: }