Actual source code: mkl_pardiso.c
petsc-3.6.4 2016-04-12
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/sbaij/seq/sbaij.h> /*I "petscmat.h" I*/
7: #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/
9: #include <stdio.h>
10: #include <stdlib.h>
11: #include <math.h>
12: #include <mkl.h>
14: /*
15: * Possible mkl_pardiso phases that controls the execution of the solver.
16: * For more information check mkl_pardiso manual.
17: */
18: #define JOB_ANALYSIS 11
19: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
20: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
21: #define JOB_NUMERICAL_FACTORIZATION 22
22: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
23: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
24: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
25: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
26: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
27: #define JOB_RELEASE_OF_LU_MEMORY 0
28: #define JOB_RELEASE_OF_ALL_MEMORY -1
30: #define IPARM_SIZE 64
32: #if defined(PETSC_USE_64BIT_INDICES)
33: #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64)
34: /* sizeof(MKL_INT) == sizeof(long long int) if ilp64*/
35: #define INT_TYPE long long int
36: #define MKL_PARDISO pardiso
37: #define MKL_PARDISO_INIT pardisoinit
38: #else
39: #define INT_TYPE long long int
40: #define MKL_PARDISO pardiso_64
41: #define MKL_PARDISO_INIT pardiso_64init
42: #endif
43: #else
44: #define INT_TYPE int
45: #define MKL_PARDISO pardiso
46: #define MKL_PARDISO_INIT pardisoinit
47: #endif
50: /*
51: * Internal data structure.
52: * For more information check mkl_pardiso manual.
53: */
54: typedef struct {
56: /* Configuration vector*/
57: INT_TYPE iparm[IPARM_SIZE];
59: /*
60: * Internal mkl_pardiso memory location.
61: * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
62: */
63: void *pt[IPARM_SIZE];
65: /* Basic mkl_pardiso info*/
66: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
68: /* Matrix structure*/
69: void *a;
70: INT_TYPE *ia, *ja;
72: /* Number of non-zero elements*/
73: INT_TYPE nz;
75: /* Row permutaton vector*/
76: INT_TYPE *perm;
78: /* Define if matrix preserves sparse structure.*/
79: MatStructure matstruc;
81: /* True if mkl_pardiso function have been used.*/
82: PetscBool CleanUp;
83: } Mat_MKL_PARDISO;
86: void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm [])
87: {
88: int iparm_copy[IPARM_SIZE], mtype_copy, i;
89:
90: mtype_copy = *mtype;
91: pardisoinit(pt, &mtype_copy, iparm_copy);
92: for(i = 0; i < IPARM_SIZE; i++){
93: iparm[i] = iparm_copy[i];
94: }
95: }
98: /*
99: * Copy the elements of matrix A.
100: * Input:
101: * - Mat A: MATSEQAIJ matrix
102: * - int shift: matrix index.
103: * - 0 for c representation
104: * - 1 for fortran representation
105: * - MatReuse reuse:
106: * - MAT_INITIAL_MATRIX: Create a new aij representation
107: * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
108: * Output:
109: * - int *nnz: Number of nonzero-elements.
110: * - int **r pointer to i index
111: * - int **c pointer to j elements
112: * - MATRIXTYPE **v: Non-zero elements
113: */
116: PetscErrorCode MatCopy_MKL_PARDISO(Mat A, MatReuse reuse, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, void **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: /*
131: * Free memory for Mat_MKL_PARDISO structure and pointers to objects.
132: */
135: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
136: {
137: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
138: PetscBool isSeqSBAIJ;
139: PetscErrorCode ierr;
142: /* Terminate instance, deallocate memories */
143: if (mat_mkl_pardiso->CleanUp) {
144: mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
146: MKL_PARDISO (mat_mkl_pardiso->pt,
147: &mat_mkl_pardiso->maxfct,
148: &mat_mkl_pardiso->mnum,
149: &mat_mkl_pardiso->mtype,
150: &mat_mkl_pardiso->phase,
151: &mat_mkl_pardiso->n,
152: NULL,
153: NULL,
154: NULL,
155: mat_mkl_pardiso->perm,
156: &mat_mkl_pardiso->nrhs,
157: mat_mkl_pardiso->iparm,
158: &mat_mkl_pardiso->msglvl,
159: NULL,
160: NULL,
161: &mat_mkl_pardiso->err);
162: }
163: PetscFree(mat_mkl_pardiso->perm);
164: PetscFree(A->spptr);
166: /* clear composed functions */
167: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
168: PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);
170: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
171: if (isSeqSBAIJ) {MatDestroy_SeqSBAIJ(A);}
172: else {MatDestroy_SeqAIJ(A);}
173: return(0);
174: }
176: /*
177: * Computes Ax = b
178: */
181: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x)
182: {
183: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
184: PetscErrorCode ierr;
185: PetscScalar *xarray;
186: const PetscScalar *barray;
189: mat_mkl_pardiso->nrhs = 1;
190: VecGetArray(x,&xarray);
191: VecGetArrayRead(b,&barray);
193: /* solve phase */
194: /*-------------*/
195: mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
196: MKL_PARDISO (mat_mkl_pardiso->pt,
197: &mat_mkl_pardiso->maxfct,
198: &mat_mkl_pardiso->mnum,
199: &mat_mkl_pardiso->mtype,
200: &mat_mkl_pardiso->phase,
201: &mat_mkl_pardiso->n,
202: mat_mkl_pardiso->a,
203: mat_mkl_pardiso->ia,
204: mat_mkl_pardiso->ja,
205: mat_mkl_pardiso->perm,
206: &mat_mkl_pardiso->nrhs,
207: mat_mkl_pardiso->iparm,
208: &mat_mkl_pardiso->msglvl,
209: (void*)barray,
210: (void*)xarray,
211: &mat_mkl_pardiso->err);
213: 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);
214: VecRestoreArray(x,&xarray);
215: VecRestoreArrayRead(b,&barray);
216: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
217: return(0);
218: }
223: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x)
224: {
225: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
226: PetscErrorCode ierr;
229: #if defined(PETSC_USE_COMPLEX)
230: mat_mkl_pardiso->iparm[12 - 1] = 1;
231: #else
232: mat_mkl_pardiso->iparm[12 - 1] = 2;
233: #endif
234: MatSolve_MKL_PARDISO(A,b,x);
235: mat_mkl_pardiso->iparm[12 - 1] = 0;
236: return(0);
237: }
242: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X)
243: {
244: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr;
245: PetscErrorCode ierr;
246: PetscScalar *barray, *xarray;
247: PetscBool flg;
250: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
251: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
252: PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
253: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");
255: MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);
257: if(mat_mkl_pardiso->nrhs > 0){
258: MatDenseGetArray(B,&barray);
259: MatDenseGetArray(X,&xarray);
261: /* solve phase */
262: /*-------------*/
263: mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
264: MKL_PARDISO (mat_mkl_pardiso->pt,
265: &mat_mkl_pardiso->maxfct,
266: &mat_mkl_pardiso->mnum,
267: &mat_mkl_pardiso->mtype,
268: &mat_mkl_pardiso->phase,
269: &mat_mkl_pardiso->n,
270: mat_mkl_pardiso->a,
271: mat_mkl_pardiso->ia,
272: mat_mkl_pardiso->ja,
273: mat_mkl_pardiso->perm,
274: &mat_mkl_pardiso->nrhs,
275: mat_mkl_pardiso->iparm,
276: &mat_mkl_pardiso->msglvl,
277: (void*)barray,
278: (void*)xarray,
279: &mat_mkl_pardiso->err);
280: 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);
281: }
282: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
283: return(0);
284: }
286: /*
287: * LU Decomposition
288: */
291: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info)
292: {
293: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->spptr;
294: PetscErrorCode ierr;
296: /* numerical factorization phase */
297: /*-------------------------------*/
299: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
300: MatCopy_MKL_PARDISO(A, MAT_REUSE_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);
302: /* numerical factorization phase */
303: /*-------------------------------*/
304: mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
305: MKL_PARDISO (mat_mkl_pardiso->pt,
306: &mat_mkl_pardiso->maxfct,
307: &mat_mkl_pardiso->mnum,
308: &mat_mkl_pardiso->mtype,
309: &mat_mkl_pardiso->phase,
310: &mat_mkl_pardiso->n,
311: mat_mkl_pardiso->a,
312: mat_mkl_pardiso->ia,
313: mat_mkl_pardiso->ja,
314: mat_mkl_pardiso->perm,
315: &mat_mkl_pardiso->nrhs,
316: mat_mkl_pardiso->iparm,
317: &mat_mkl_pardiso->msglvl,
318: NULL,
319: NULL,
320: &mat_mkl_pardiso->err);
321: 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);
323: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
324: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
325: return(0);
326: }
328: /* Sets mkl_pardiso options from the options database */
331: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A)
332: {
333: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
334: PetscErrorCode ierr;
335: PetscInt icntl, threads = 1;
336: PetscBool flg;
339: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");
340: PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use","None",threads,&threads,&flg);
341: if (flg) mkl_set_num_threads((int)threads);
343: 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);
344: if (flg) mat_mkl_pardiso->maxfct = icntl;
346: PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);
347: if (flg) mat_mkl_pardiso->mnum = icntl;
348:
349: PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);
350: if (flg) mat_mkl_pardiso->msglvl = icntl;
352: PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);
353: if(flg){
354: void *pt[IPARM_SIZE];
355: mat_mkl_pardiso->mtype = icntl;
356: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
357: #if defined(PETSC_USE_REAL_SINGLE)
358: mat_mkl_pardiso->iparm[27] = 1;
359: #else
360: mat_mkl_pardiso->iparm[27] = 0;
361: #endif
362: mat_mkl_pardiso->iparm[34] = 1;
363: }
364: PetscOptionsInt("-mat_mkl_pardiso_1","Use default values","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);
366: if(flg && icntl != 0){
367: PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);
368: if (flg) mat_mkl_pardiso->iparm[1] = icntl;
370: PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);
371: if (flg) mat_mkl_pardiso->iparm[3] = icntl;
373: PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);
374: if (flg) mat_mkl_pardiso->iparm[4] = icntl;
376: PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);
377: if (flg) mat_mkl_pardiso->iparm[5] = icntl;
379: PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);
380: if (flg) mat_mkl_pardiso->iparm[7] = icntl;
382: PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);
383: if (flg) mat_mkl_pardiso->iparm[9] = icntl;
385: PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);
386: if (flg) mat_mkl_pardiso->iparm[10] = icntl;
388: PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);
389: if (flg) mat_mkl_pardiso->iparm[11] = icntl;
391: PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);
392: if (flg) mat_mkl_pardiso->iparm[12] = icntl;
394: PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);
395: if (flg) mat_mkl_pardiso->iparm[17] = icntl;
397: PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);
398: if (flg) mat_mkl_pardiso->iparm[18] = icntl;
400: PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);
401: if (flg) mat_mkl_pardiso->iparm[20] = icntl;
403: PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);
404: if (flg) mat_mkl_pardiso->iparm[23] = icntl;
406: PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);
407: if (flg) mat_mkl_pardiso->iparm[24] = icntl;
409: PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);
410: if (flg) mat_mkl_pardiso->iparm[26] = icntl;
412: PetscOptionsInt("-mat_mkl_pardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_pardiso->iparm[30],&icntl,&flg);
413: if (flg) mat_mkl_pardiso->iparm[30] = icntl;
415: PetscOptionsInt("-mat_mkl_pardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_pardiso->iparm[33],&icntl,&flg);
416: if (flg) mat_mkl_pardiso->iparm[33] = icntl;
418: PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
419: if (flg) mat_mkl_pardiso->iparm[59] = icntl;
420: }
421: PetscOptionsEnd();
422: return(0);
423: }
427: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
428: {
430: PetscInt i;
433: for ( i = 0; i < IPARM_SIZE; i++ ){
434: mat_mkl_pardiso->iparm[i] = 0;
435: }
437: for ( i = 0; i < IPARM_SIZE; i++ ){
438: mat_mkl_pardiso->pt[i] = 0;
439: }
440:
441: /* Default options for both sym and unsym */
442: mat_mkl_pardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */
443: mat_mkl_pardiso->iparm[ 1] = 2; /* Metis reordering */
444: mat_mkl_pardiso->iparm[ 5] = 0; /* Write solution into x */
445: mat_mkl_pardiso->iparm[ 7] = 2; /* Max number of iterative refinement steps */
446: mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
447: mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
448: #if 0
449: mat_mkl_pardiso->iparm[23] = 1; /* Parallel factorization control*/
450: #endif
451: mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
452: mat_mkl_pardiso->iparm[39] = 0; /* Input: matrix/rhs/solution stored on master */
453:
454: mat_mkl_pardiso->CleanUp = PETSC_FALSE;
455: mat_mkl_pardiso->maxfct = 1; /* Maximum number of numerical factorizations. */
456: mat_mkl_pardiso->mnum = 1; /* Which factorization to use. */
457: mat_mkl_pardiso->msglvl = 0; /* 0: do not print 1: Print statistical information in file */
458: mat_mkl_pardiso->phase = -1;
459: mat_mkl_pardiso->err = 0;
460:
461: mat_mkl_pardiso->n = A->rmap->N;
462: mat_mkl_pardiso->nrhs = 1;
463: mat_mkl_pardiso->err = 0;
464: mat_mkl_pardiso->phase = -1;
465:
466: if(ftype == MAT_FACTOR_LU){
467: /* Default type for non-sym */
468: #if defined(PETSC_USE_COMPLEX)
469: mat_mkl_pardiso->mtype = 13;
470: #else
471: mat_mkl_pardiso->mtype = 11;
472: #endif
474: mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
475: mat_mkl_pardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
476: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
478: } else {
479: /* Default type for sym */
480: #if defined(PETSC_USE_COMPLEX)
481: mat_mkl_pardiso ->mtype = 3;
482: #else
483: mat_mkl_pardiso ->mtype = -2;
484: #endif
485: mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
486: mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
487: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
488: /* mat_mkl_pardiso->iparm[20] = 1; */ /* Apply 1x1 and 2x2 Bunch-Kaufman pivoting during the factorization process */
489: #if defined(PETSC_USE_DEBUG)
490: mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
491: #endif
492: }
493: PetscMalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
494: for(i = 0; i < A->rmap->N; i++){
495: mat_mkl_pardiso->perm[i] = 0;
496: }
497: return(0);
498: }
500: /*
501: * Symbolic decomposition. Mkl_Pardiso analysis phase.
502: */
505: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
506: {
507: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr;
508: PetscErrorCode ierr;
511: mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
513: /* Set MKL_PARDISO options from the options database */
514: PetscSetMKL_PARDISOFromOptions(F,A);
516: MatCopy_MKL_PARDISO(A, MAT_INITIAL_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);
517: mat_mkl_pardiso->n = A->rmap->N;
519: /* analysis phase */
520: /*----------------*/
521: mat_mkl_pardiso->phase = JOB_ANALYSIS;
523: MKL_PARDISO (mat_mkl_pardiso->pt,
524: &mat_mkl_pardiso->maxfct,
525: &mat_mkl_pardiso->mnum,
526: &mat_mkl_pardiso->mtype,
527: &mat_mkl_pardiso->phase,
528: &mat_mkl_pardiso->n,
529: mat_mkl_pardiso->a,
530: mat_mkl_pardiso->ia,
531: mat_mkl_pardiso->ja,
532: mat_mkl_pardiso->perm,
533: &mat_mkl_pardiso->nrhs,
534: mat_mkl_pardiso->iparm,
535: &mat_mkl_pardiso->msglvl,
536: NULL,
537: NULL,
538: &mat_mkl_pardiso->err);
539: 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);
541: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
543: if(F->factortype == MAT_FACTOR_LU){
544: F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
545: } else {
546: F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO;
547: }
548: F->ops->solve = MatSolve_MKL_PARDISO;
549: F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO;
550: F->ops->matsolve = MatMatSolve_MKL_PARDISO;
551: return(0);
552: }
556: PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
557: {
561: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
562: return(0);
563: }
567: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,const MatFactorInfo *info)
568: {
572: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
573: return(0);
574: }
578: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
579: {
580: PetscErrorCode ierr;
581: PetscBool iascii;
582: PetscViewerFormat format;
583: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr;
584: PetscInt i;
587: /* check if matrix is mkl_pardiso type */
588: if (A->ops->solve != MatSolve_MKL_PARDISO) return(0);
590: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
591: if (iascii) {
592: PetscViewerGetFormat(viewer,&format);
593: if (format == PETSC_VIEWER_ASCII_INFO) {
594: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
595: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);
596: for(i = 1; i <= 64; i++){
597: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
598: }
599: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);
600: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);
601: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);
602: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);
603: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);
604: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);
605: }
606: }
607: return(0);
608: }
613: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
614: {
615: Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->spptr;
618: info->block_size = 1.0;
619: info->nz_allocated = mat_mkl_pardiso->nz + 0.0;
620: info->nz_unneeded = 0.0;
621: info->assemblies = 0.0;
622: info->mallocs = 0.0;
623: info->memory = 0.0;
624: info->fill_ratio_given = 0;
625: info->fill_ratio_needed = 0;
626: info->factor_mallocs = 0;
627: return(0);
628: }
632: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
633: {
634: Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->spptr;
637: if(icntl <= 64){
638: mat_mkl_pardiso->iparm[icntl - 1] = ival;
639: } else {
640: if(icntl == 65)
641: mkl_set_num_threads((int)ival);
642: else if(icntl == 66)
643: mat_mkl_pardiso->maxfct = ival;
644: else if(icntl == 67)
645: mat_mkl_pardiso->mnum = ival;
646: else if(icntl == 68)
647: mat_mkl_pardiso->msglvl = ival;
648: else if(icntl == 69){
649: void *pt[IPARM_SIZE];
650: mat_mkl_pardiso->mtype = ival;
651: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
652: #if defined(PETSC_USE_REAL_SINGLE)
653: mat_mkl_pardiso->iparm[27] = 1;
654: #else
655: mat_mkl_pardiso->iparm[27] = 0;
656: #endif
657: mat_mkl_pardiso->iparm[34] = 1;
658: }
659: }
660: return(0);
661: }
665: /*@
666: MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters
668: Logically Collective on Mat
670: Input Parameters:
671: + F - the factored matrix obtained by calling MatGetFactor()
672: . icntl - index of Mkl_Pardiso parameter
673: - ival - value of Mkl_Pardiso parameter
675: Options Database:
676: . -mat_mkl_pardiso_<icntl> <ival>
678: Level: beginner
680: References: Mkl_Pardiso Users' Guide
682: .seealso: MatGetFactor()
683: @*/
684: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
685: {
689: PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
690: return(0);
691: }
693: /*MC
694: MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for
695: sequential matrices via the external package MKL_PARDISO.
697: Works with MATSEQAIJ matrices
699: Use -pc_type lu -pc_factor_mat_solver_package mkl_pardiso to us this direct solver
701: Options Database Keys:
702: + -mat_mkl_pardiso_65 - Number of threads to use
703: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
704: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
705: . -mat_mkl_pardiso_68 - Message level information
706: . -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
707: . -mat_mkl_pardiso_1 - Use default values
708: . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix
709: . -mat_mkl_pardiso_4 - Preconditioned CGS/CG
710: . -mat_mkl_pardiso_5 - User permutation
711: . -mat_mkl_pardiso_6 - Write solution on x
712: . -mat_mkl_pardiso_8 - Iterative refinement step
713: . -mat_mkl_pardiso_10 - Pivoting perturbation
714: . -mat_mkl_pardiso_11 - Scaling vectors
715: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
716: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
717: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
718: . -mat_mkl_pardiso_19 - Report number of floating point operations
719: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
720: . -mat_mkl_pardiso_24 - Parallel factorization control
721: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
722: . -mat_mkl_pardiso_27 - Matrix checker
723: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
724: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
725: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode
727: Level: beginner
729: For more information please check mkl_pardiso manual
731: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
733: M*/
736: static PetscErrorCode MatFactorGetSolverPackage_mkl_pardiso(Mat A, const MatSolverPackage *type)
737: {
739: *type = MATSOLVERMKL_PARDISO;
740: return(0);
741: }
743: /* MatGetFactor for Seq sbAIJ matrices */
746: PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
747: {
748: Mat B;
749: PetscErrorCode ierr;
750: Mat_MKL_PARDISO *mat_mkl_pardiso;
751: PetscBool isSeqSBAIJ;
752: PetscInt bs;
755: /* Create the factorization matrix */
756: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
757: MatCreate(PetscObjectComm((PetscObject)A),&B);
758: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
759: MatSetType(B,((PetscObject)A)->type_name);
760: MatSeqSBAIJSetPreallocation(B,1,0,NULL);
761: MatGetBlockSize(A,&bs);
763: if(bs != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQSBAIJ with block size other than 1 is not supported by Pardiso");
764: if(ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_CHOLESKY.");
765:
766: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
767: B->factortype = MAT_FACTOR_CHOLESKY;
768: B->ops->destroy = MatDestroy_MKL_PARDISO;
769: B->ops->view = MatView_MKL_PARDISO;
770: B->factortype = ftype;
771: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
772: B->assembled = PETSC_TRUE; /* required by -ksp_view */
774: PetscNewLog(B,&mat_mkl_pardiso);
775: B->spptr = mat_mkl_pardiso;
776: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);
777: PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
778: MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
779: *F = B;
780: return(0);
781: }
783: /* MatGetFactor for Seq AIJ matrices */
786: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
787: {
788: Mat B;
789: PetscErrorCode ierr;
790: Mat_MKL_PARDISO *mat_mkl_pardiso;
791: PetscBool isSeqAIJ;
794: /* Create the factorization matrix */
795: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
796: MatCreate(PetscObjectComm((PetscObject)A),&B);
797: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
798: MatSetType(B,((PetscObject)A)->type_name);
799: MatSeqAIJSetPreallocation(B,0,NULL);
801: if(ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_LU.");
803: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
804: B->factortype = MAT_FACTOR_LU;
805: B->ops->destroy = MatDestroy_MKL_PARDISO;
806: B->ops->view = MatView_MKL_PARDISO;
807: B->factortype = ftype;
808: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
809: B->assembled = PETSC_TRUE; /* required by -ksp_view */
811: PetscNewLog(B,&mat_mkl_pardiso);
812: B->spptr = mat_mkl_pardiso;
813: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);
814: PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
815: MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
817: *F = B;
818: return(0);
819: }
823: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MKL_Pardiso(void)
824: {
828: MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_mkl_pardiso );
829: MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mkl_pardiso);
830: return(0);
831: }