Actual source code: mkl_pardiso.c
petsc-3.6.1 2015-08-06
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/dense/seq/dense.h>
8: #include <stdio.h>
9: #include <stdlib.h>
10: #include <math.h>
11: #include <mkl.h>
13: /*
14: * Possible mkl_pardiso phases that controls the execution of the solver.
15: * For more information check mkl_pardiso manual.
16: */
17: #define JOB_ANALYSIS 11
18: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
19: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
20: #define JOB_NUMERICAL_FACTORIZATION 22
21: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
22: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
23: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
24: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
25: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
26: #define JOB_RELEASE_OF_LU_MEMORY 0
27: #define JOB_RELEASE_OF_ALL_MEMORY -1
29: #define IPARM_SIZE 64
31: #if defined(PETSC_USE_64BIT_INDICES)
32: #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64)
33: /* sizeof(MKL_INT) == sizeof(long long int) if ilp64*/
34: #define INT_TYPE long long int
35: #define MKL_PARDISO pardiso
36: #define MKL_PARDISO_INIT pardisoinit
37: #else
38: #define INT_TYPE long long int
39: #define MKL_PARDISO pardiso_64
40: #define MKL_PARDISO_INIT pardiso_64init
41: #endif
42: #else
43: #define INT_TYPE int
44: #define MKL_PARDISO pardiso
45: #define MKL_PARDISO_INIT pardisoinit
46: #endif
49: /*
50: * Internal data structure.
51: * For more information check mkl_pardiso manual.
52: */
53: typedef struct {
55: /* Configuration vector*/
56: INT_TYPE iparm[IPARM_SIZE];
58: /*
59: * Internal mkl_pardiso memory location.
60: * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
61: */
62: void *pt[IPARM_SIZE];
64: /* Basic mkl_pardiso info*/
65: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
67: /* Matrix structure*/
68: void *a;
69: INT_TYPE *ia, *ja;
71: /* Number of non-zero elements*/
72: INT_TYPE nz;
74: /* Row permutaton vector*/
75: INT_TYPE *perm;
77: /* Define if matrix preserves sparse structure.*/
78: MatStructure matstruc;
80: /* True if mkl_pardiso function have been used.*/
81: PetscBool CleanUp;
82: } Mat_MKL_PARDISO;
85: void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm [])
86: {
87: int iparm_copy[IPARM_SIZE], mtype_copy, i;
88:
89: mtype_copy = *mtype;
90: pardisoinit(pt, &mtype_copy, iparm_copy);
91: for(i = 0; i < IPARM_SIZE; i++){
92: iparm[i] = iparm_copy[i];
93: }
94: }
97: /*
98: * Copy the elements of matrix A.
99: * Input:
100: * - Mat A: MATSEQAIJ matrix
101: * - int shift: matrix index.
102: * - 0 for c representation
103: * - 1 for fortran representation
104: * - MatReuse reuse:
105: * - MAT_INITIAL_MATRIX: Create a new aij representation
106: * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
107: * Output:
108: * - int *nnz: Number of nonzero-elements.
109: * - int **r pointer to i index
110: * - int **c pointer to j elements
111: * - MATRIXTYPE **v: Non-zero elements
112: */
115: PetscErrorCode MatCopy_MKL_PARDISO(Mat A, MatReuse reuse, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, void **v)
116: {
117: Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data;
120: *v=aa->a;
121: if (reuse == MAT_INITIAL_MATRIX) {
122: *r = (INT_TYPE*)aa->i;
123: *c = (INT_TYPE*)aa->j;
124: *nnz = aa->nz;
125: }
126: return(0);
127: }
129: /*
130: * Free memory for Mat_MKL_PARDISO structure and pointers to objects.
131: */
134: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
135: {
136: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
137: PetscErrorCode ierr;
140: /* Terminate instance, deallocate memories */
141: if (mat_mkl_pardiso->CleanUp) {
142: mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
144: MKL_PARDISO (mat_mkl_pardiso->pt,
145: &mat_mkl_pardiso->maxfct,
146: &mat_mkl_pardiso->mnum,
147: &mat_mkl_pardiso->mtype,
148: &mat_mkl_pardiso->phase,
149: &mat_mkl_pardiso->n,
150: NULL,
151: NULL,
152: NULL,
153: mat_mkl_pardiso->perm,
154: &mat_mkl_pardiso->nrhs,
155: mat_mkl_pardiso->iparm,
156: &mat_mkl_pardiso->msglvl,
157: NULL,
158: NULL,
159: &mat_mkl_pardiso->err);
160: }
161: PetscFree(mat_mkl_pardiso->perm);
162: PetscFree(A->spptr);
164: /* clear composed functions */
165: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
166: PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);
168: MatDestroy_SeqAIJ(A);
169: return(0);
170: }
172: /*
173: * Computes Ax = b
174: */
177: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x)
178: {
179: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
180: PetscErrorCode ierr;
181: PetscScalar *xarray;
182: const PetscScalar *barray;
185: mat_mkl_pardiso->nrhs = 1;
186: VecGetArray(x,&xarray);
187: VecGetArrayRead(b,&barray);
189: /* solve phase */
190: /*-------------*/
191: mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
192: MKL_PARDISO (mat_mkl_pardiso->pt,
193: &mat_mkl_pardiso->maxfct,
194: &mat_mkl_pardiso->mnum,
195: &mat_mkl_pardiso->mtype,
196: &mat_mkl_pardiso->phase,
197: &mat_mkl_pardiso->n,
198: mat_mkl_pardiso->a,
199: mat_mkl_pardiso->ia,
200: mat_mkl_pardiso->ja,
201: mat_mkl_pardiso->perm,
202: &mat_mkl_pardiso->nrhs,
203: mat_mkl_pardiso->iparm,
204: &mat_mkl_pardiso->msglvl,
205: (void*)barray,
206: (void*)xarray,
207: &mat_mkl_pardiso->err);
209: 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);
210: VecRestoreArray(x,&xarray);
211: VecRestoreArrayRead(b,&barray);
212: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
213: return(0);
214: }
219: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x)
220: {
221: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
222: PetscErrorCode ierr;
225: #if defined(PETSC_USE_COMPLEX)
226: mat_mkl_pardiso->iparm[12 - 1] = 1;
227: #else
228: mat_mkl_pardiso->iparm[12 - 1] = 2;
229: #endif
230: MatSolve_MKL_PARDISO(A,b,x);
231: mat_mkl_pardiso->iparm[12 - 1] = 0;
232: return(0);
233: }
238: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X)
239: {
240: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
241: PetscErrorCode ierr;
242: PetscScalar *barray, *xarray;
243: PetscBool flg;
246: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
247: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
248: PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
249: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");
251: MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);
253: if(mat_mkl_pardiso->nrhs > 0){
254: MatDenseGetArray(B,&barray);
255: MatDenseGetArray(X,&xarray);
257: /* solve phase */
258: /*-------------*/
259: mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
260: MKL_PARDISO (mat_mkl_pardiso->pt,
261: &mat_mkl_pardiso->maxfct,
262: &mat_mkl_pardiso->mnum,
263: &mat_mkl_pardiso->mtype,
264: &mat_mkl_pardiso->phase,
265: &mat_mkl_pardiso->n,
266: mat_mkl_pardiso->a,
267: mat_mkl_pardiso->ia,
268: mat_mkl_pardiso->ja,
269: mat_mkl_pardiso->perm,
270: &mat_mkl_pardiso->nrhs,
271: mat_mkl_pardiso->iparm,
272: &mat_mkl_pardiso->msglvl,
273: (void*)barray,
274: (void*)xarray,
275: &mat_mkl_pardiso->err);
276: 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);
277: }
278: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
279: return(0);
280: }
282: /*
283: * LU Decomposition
284: */
287: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info)
288: {
289: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->spptr;
290: PetscErrorCode ierr;
292: /* numerical factorization phase */
293: /*-------------------------------*/
295: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
296: MatCopy_MKL_PARDISO(A, MAT_REUSE_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);
298: /* numerical factorization phase */
299: /*-------------------------------*/
300: mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
301: MKL_PARDISO (mat_mkl_pardiso->pt,
302: &mat_mkl_pardiso->maxfct,
303: &mat_mkl_pardiso->mnum,
304: &mat_mkl_pardiso->mtype,
305: &mat_mkl_pardiso->phase,
306: &mat_mkl_pardiso->n,
307: mat_mkl_pardiso->a,
308: mat_mkl_pardiso->ia,
309: mat_mkl_pardiso->ja,
310: mat_mkl_pardiso->perm,
311: &mat_mkl_pardiso->nrhs,
312: mat_mkl_pardiso->iparm,
313: &mat_mkl_pardiso->msglvl,
314: NULL,
315: NULL,
316: &mat_mkl_pardiso->err);
317: 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);
319: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
320: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
321: return(0);
322: }
324: /* Sets mkl_pardiso options from the options database */
327: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A)
328: {
329: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
330: PetscErrorCode ierr;
331: PetscInt icntl;
332: PetscBool flg;
333: int pt[IPARM_SIZE], threads = 1;
336: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");
337: PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use","None",threads,&threads,&flg);
338: if (flg) mkl_set_num_threads(threads);
340: 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);
341: if (flg) mat_mkl_pardiso->maxfct = icntl;
343: PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);
344: if (flg) mat_mkl_pardiso->mnum = icntl;
345:
346: PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);
347: if (flg) mat_mkl_pardiso->msglvl = icntl;
349: PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);
350: if(flg){
351: mat_mkl_pardiso->mtype = icntl;
352: MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
353: #if defined(PETSC_USE_REAL_SINGLE)
354: mat_mkl_pardiso->iparm[27] = 1;
355: #else
356: mat_mkl_pardiso->iparm[27] = 0;
357: #endif
358: mat_mkl_pardiso->iparm[34] = 1;
359: }
360: PetscOptionsInt("-mat_mkl_pardiso_1","Use default values","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);
362: if(flg && icntl != 0){
363: PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);
364: if (flg) mat_mkl_pardiso->iparm[1] = icntl;
366: PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);
367: if (flg) mat_mkl_pardiso->iparm[3] = icntl;
369: PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);
370: if (flg) mat_mkl_pardiso->iparm[4] = icntl;
372: PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);
373: if (flg) mat_mkl_pardiso->iparm[5] = icntl;
375: PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);
376: if (flg) mat_mkl_pardiso->iparm[7] = icntl;
378: PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);
379: if (flg) mat_mkl_pardiso->iparm[9] = icntl;
381: PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);
382: if (flg) mat_mkl_pardiso->iparm[10] = icntl;
384: PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);
385: if (flg) mat_mkl_pardiso->iparm[11] = icntl;
387: PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);
388: if (flg) mat_mkl_pardiso->iparm[12] = icntl;
390: PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);
391: if (flg) mat_mkl_pardiso->iparm[17] = icntl;
393: PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);
394: if (flg) mat_mkl_pardiso->iparm[18] = icntl;
396: PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);
397: if (flg) mat_mkl_pardiso->iparm[20] = icntl;
399: PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);
400: if (flg) mat_mkl_pardiso->iparm[23] = icntl;
402: PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);
403: if (flg) mat_mkl_pardiso->iparm[24] = icntl;
405: PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);
406: if (flg) mat_mkl_pardiso->iparm[26] = icntl;
408: PetscOptionsInt("-mat_mkl_pardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_pardiso->iparm[30],&icntl,&flg);
409: if (flg) mat_mkl_pardiso->iparm[30] = icntl;
411: PetscOptionsInt("-mat_mkl_pardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_pardiso->iparm[33],&icntl,&flg);
412: if (flg) mat_mkl_pardiso->iparm[33] = icntl;
414: PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
415: if (flg) mat_mkl_pardiso->iparm[59] = icntl;
416: }
417: PetscOptionsEnd();
418: return(0);
419: }
423: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
424: {
426: PetscInt i;
429: for ( i = 0; i < IPARM_SIZE; i++ ){
430: mat_mkl_pardiso->iparm[i] = 0;
431: }
433: for ( i = 0; i < IPARM_SIZE; i++ ){
434: mat_mkl_pardiso->pt[i] = 0;
435: }
436:
437: /* Default options for both sym and unsym */
438: mat_mkl_pardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */
439: mat_mkl_pardiso->iparm[ 1] = 2; /* Metis reordering */
440: mat_mkl_pardiso->iparm[ 5] = 0; /* Write solution into x */
441: mat_mkl_pardiso->iparm[ 7] = 2; /* Max number of iterative refinement steps */
442: mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
443: mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
444: #if 0
445: mat_mkl_pardiso->iparm[23] = 1; /* Parallel factorization control*/
446: #endif
447: mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
448: mat_mkl_pardiso->iparm[39] = 0; /* Input: matrix/rhs/solution stored on master */
449:
450: mat_mkl_pardiso->CleanUp = PETSC_FALSE;
451: mat_mkl_pardiso->maxfct = 1; /* Maximum number of numerical factorizations. */
452: mat_mkl_pardiso->mnum = 1; /* Which factorization to use. */
453: mat_mkl_pardiso->msglvl = 0; /* 0: do not print 1: Print statistical information in file */
454: mat_mkl_pardiso->phase = -1;
455: mat_mkl_pardiso->err = 0;
456:
457: mat_mkl_pardiso->n = A->rmap->N;
458: mat_mkl_pardiso->nrhs = 1;
459: mat_mkl_pardiso->err = 0;
460: mat_mkl_pardiso->phase = -1;
461:
462: if(ftype == MAT_FACTOR_LU){
463: /* Default type for non-sym */
464: #if defined(PETSC_USE_COMPLEX)
465: mat_mkl_pardiso->mtype = 13;
466: #else
467: mat_mkl_pardiso->mtype = 11;
468: #endif
470: mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
471: mat_mkl_pardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
472: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
474: } else {
475: /* Default type for sym */
476: #if defined(PETSC_USE_COMPLEX)
477: mat_mkl_pardiso ->mtype = 3;
478: #else
479: mat_mkl_pardiso ->mtype = -2;
480: #endif
481: mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
482: mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
483: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
484: /* mat_mkl_pardiso->iparm[20] = 1; */ /* Apply 1x1 and 2x2 Bunch-Kaufman pivoting during the factorization process */
485: #if defined(PETSC_USE_DEBUG)
486: mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
487: #endif
488: }
489: PetscMalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
490: for(i = 0; i < A->rmap->N; i++){
491: mat_mkl_pardiso->perm[i] = 0;
492: }
493: return(0);
494: }
496: /*
497: * Symbolic decomposition. Mkl_Pardiso analysis phase.
498: */
501: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
502: {
503: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
504: PetscErrorCode ierr;
507: mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
509: /* Set MKL_PARDISO options from the options database */
510: PetscSetMKL_PARDISOFromOptions(F,A);
512: MatCopy_MKL_PARDISO(A, MAT_INITIAL_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);
513: mat_mkl_pardiso->n = A->rmap->N;
515: /* analysis phase */
516: /*----------------*/
517: mat_mkl_pardiso->phase = JOB_ANALYSIS;
519: MKL_PARDISO (mat_mkl_pardiso->pt,
520: &mat_mkl_pardiso->maxfct,
521: &mat_mkl_pardiso->mnum,
522: &mat_mkl_pardiso->mtype,
523: &mat_mkl_pardiso->phase,
524: &mat_mkl_pardiso->n,
525: mat_mkl_pardiso->a,
526: mat_mkl_pardiso->ia,
527: mat_mkl_pardiso->ja,
528: mat_mkl_pardiso->perm,
529: &mat_mkl_pardiso->nrhs,
530: mat_mkl_pardiso->iparm,
531: &mat_mkl_pardiso->msglvl,
532: NULL,
533: NULL,
534: &mat_mkl_pardiso->err);
535: 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);
537: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
539: if(F->factortype == MAT_FACTOR_LU){
540: F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
541: } else {
542: F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO;
543: }
544: F->ops->solve = MatSolve_MKL_PARDISO;
545: F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO;
546: F->ops->matsolve = MatMatSolve_MKL_PARDISO;
547: return(0);
548: }
552: PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
553: {
557: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
558: return(0);
559: }
563: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,const MatFactorInfo *info)
564: {
568: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
569: return(0);
570: }
574: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
575: {
576: PetscErrorCode ierr;
577: PetscBool iascii;
578: PetscViewerFormat format;
579: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
580: PetscInt i;
583: /* check if matrix is mkl_pardiso type */
584: if (A->ops->solve != MatSolve_MKL_PARDISO) return(0);
586: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
587: if (iascii) {
588: PetscViewerGetFormat(viewer,&format);
589: if (format == PETSC_VIEWER_ASCII_INFO) {
590: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
591: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);
592: for(i = 1; i <= 64; i++){
593: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
594: }
595: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);
596: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);
597: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);
598: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);
599: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);
600: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);
601: }
602: }
603: return(0);
604: }
609: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
610: {
611: Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->spptr;
614: info->block_size = 1.0;
615: info->nz_allocated = mat_mkl_pardiso->nz + 0.0;
616: info->nz_unneeded = 0.0;
617: info->assemblies = 0.0;
618: info->mallocs = 0.0;
619: info->memory = 0.0;
620: info->fill_ratio_given = 0;
621: info->fill_ratio_needed = 0;
622: info->factor_mallocs = 0;
623: return(0);
624: }
628: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
629: {
630: Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->spptr;
633: if(icntl <= 64){
634: mat_mkl_pardiso->iparm[icntl - 1] = ival;
635: } else {
636: if(icntl == 65)
637: mkl_set_num_threads((int)ival);
638: else if(icntl == 66)
639: mat_mkl_pardiso->maxfct = ival;
640: else if(icntl == 67)
641: mat_mkl_pardiso->mnum = ival;
642: else if(icntl == 68)
643: mat_mkl_pardiso->msglvl = ival;
644: else if(icntl == 69){
645: int pt[IPARM_SIZE];
646: mat_mkl_pardiso->mtype = ival;
647: MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
648: #if defined(PETSC_USE_REAL_SINGLE)
649: mat_mkl_pardiso->iparm[27] = 1;
650: #else
651: mat_mkl_pardiso->iparm[27] = 0;
652: #endif
653: mat_mkl_pardiso->iparm[34] = 1;
654: }
655: }
656: return(0);
657: }
661: /*@
662: MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters
664: Logically Collective on Mat
666: Input Parameters:
667: + F - the factored matrix obtained by calling MatGetFactor()
668: . icntl - index of Mkl_Pardiso parameter
669: - ival - value of Mkl_Pardiso parameter
671: Options Database:
672: . -mat_mkl_pardiso_<icntl> <ival>
674: Level: beginner
676: References: Mkl_Pardiso Users' Guide
678: .seealso: MatGetFactor()
679: @*/
680: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
681: {
685: PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
686: return(0);
687: }
689: /*MC
690: MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for
691: sequential matrices via the external package MKL_PARDISO.
693: Works with MATSEQAIJ matrices
695: Use -pc_type lu -pc_factor_mat_solver_package mkl_pardiso to us this direct solver
697: Options Database Keys:
698: + -mat_mkl_pardiso_65 - Number of thrads to use
699: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
700: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
701: . -mat_mkl_pardiso_68 - Message level information
702: . -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
703: . -mat_mkl_pardiso_1 - Use default values
704: . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix
705: . -mat_mkl_pardiso_4 - Preconditioned CGS/CG
706: . -mat_mkl_pardiso_5 - User permutation
707: . -mat_mkl_pardiso_6 - Write solution on x
708: . -mat_mkl_pardiso_8 - Iterative refinement step
709: . -mat_mkl_pardiso_10 - Pivoting perturbation
710: . -mat_mkl_pardiso_11 - Scaling vectors
711: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
712: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
713: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
714: . -mat_mkl_pardiso_19 - Report number of floating point operations
715: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
716: . -mat_mkl_pardiso_24 - Parallel factorization control
717: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
718: . -mat_mkl_pardiso_27 - Matrix checker
719: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
720: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
721: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode
723: Level: beginner
725: For more information please check mkl_pardiso manual
727: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
729: M*/
732: static PetscErrorCode MatFactorGetSolverPackage_mkl_pardiso(Mat A, const MatSolverPackage *type)
733: {
735: *type = MATSOLVERMKL_PARDISO;
736: return(0);
737: }
739: /* MatGetFactor for Seq sbAIJ matrices */
742: PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
743: {
744: Mat B;
745: PetscErrorCode ierr;
746: Mat_MKL_PARDISO *mat_mkl_pardiso;
747: PetscBool isSeqSBAIJ;
748: PetscInt bs;
751: /* Create the factorization matrix */
752: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
753: MatCreate(PetscObjectComm((PetscObject)A),&B);
754: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
755: MatSetType(B,((PetscObject)A)->type_name);
756: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
757: MatGetBlockSize(A,&bs);
759: if(bs != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQSBAIJ with block size other than 1 is not supported by Pardiso");
760: if(ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_CHOLESKY.");
761:
762: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
763: B->factortype = MAT_FACTOR_CHOLESKY;
764: B->ops->destroy = MatDestroy_MKL_PARDISO;
765: B->ops->view = MatView_MKL_PARDISO;
766: B->factortype = ftype;
767: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
768: B->assembled = PETSC_TRUE; /* required by -ksp_view */
770: PetscNewLog(B,&mat_mkl_pardiso);
771: B->spptr = mat_mkl_pardiso;
772: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);
773: PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
774: MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
775: *F = B;
776: return(0);
777: }
779: /* MatGetFactor for Seq AIJ matrices */
782: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
783: {
784: Mat B;
785: PetscErrorCode ierr;
786: Mat_MKL_PARDISO *mat_mkl_pardiso;
787: PetscBool isSeqAIJ;
790: /* Create the factorization matrix */
791: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
792: MatCreate(PetscObjectComm((PetscObject)A),&B);
793: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
794: MatSetType(B,((PetscObject)A)->type_name);
795: MatSeqAIJSetPreallocation(B,0,NULL);
797: if(ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_LU.");
799: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
800: B->factortype = MAT_FACTOR_LU;
801: B->ops->destroy = MatDestroy_MKL_PARDISO;
802: B->ops->view = MatView_MKL_PARDISO;
803: B->factortype = ftype;
804: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
805: B->assembled = PETSC_TRUE; /* required by -ksp_view */
807: PetscNewLog(B,&mat_mkl_pardiso);
808: B->spptr = mat_mkl_pardiso;
809: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);
810: PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
811: MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
813: *F = B;
814: return(0);
815: }
819: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MKL_Pardiso(void)
820: {
824: MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_mkl_pardiso );
825: MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mkl_pardiso);
826: return(0);
827: }