Actual source code: mkl_cpardiso.c
petsc-3.12.5 2020-03-29
2: #include <petscsys.h>
3: #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/
4: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
6: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
7: #define MKL_ILP64
8: #endif
9: #include <mkl.h>
10: #include <mkl_cluster_sparse_solver.h>
12: /*
13: * Possible mkl_cpardiso phases that controls the execution of the solver.
14: * For more information check mkl_cpardiso manual.
15: */
16: #define JOB_ANALYSIS 11
17: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
18: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
19: #define JOB_NUMERICAL_FACTORIZATION 22
20: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
21: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
22: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
23: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
24: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
25: #define JOB_RELEASE_OF_LU_MEMORY 0
26: #define JOB_RELEASE_OF_ALL_MEMORY -1
28: #define IPARM_SIZE 64
29: #define INT_TYPE MKL_INT
31: static const char *Err_MSG_CPardiso(int errNo) {
32: switch (errNo) {
33: case -1:
34: return "input inconsistent"; break;
35: case -2:
36: return "not enough memory"; break;
37: case -3:
38: return "reordering problem"; break;
39: case -4:
40: return "zero pivot, numerical factorization or iterative refinement problem"; break;
41: case -5:
42: return "unclassified (internal) error"; break;
43: case -6:
44: return "preordering failed (matrix types 11, 13 only)"; break;
45: case -7:
46: return "diagonal matrix problem"; break;
47: case -8:
48: return "32-bit integer overflow problem"; break;
49: case -9:
50: return "not enough memory for OOC"; break;
51: case -10:
52: return "problems with opening OOC temporary files"; break;
53: case -11:
54: return "read/write problems with the OOC data file"; break;
55: default :
56: return "unknown error";
57: }
58: }
60: /*
61: * Internal data structure.
62: * For more information check mkl_cpardiso manual.
63: */
65: typedef struct {
67: /* Configuration vector */
68: INT_TYPE iparm[IPARM_SIZE];
70: /*
71: * Internal mkl_cpardiso memory location.
72: * After the first call to mkl_cpardiso do not modify pt, as that could cause a serious memory leak.
73: */
74: void *pt[IPARM_SIZE];
76: MPI_Comm comm_mkl_cpardiso;
78: /* Basic mkl_cpardiso info*/
79: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
81: /* Matrix structure */
82: PetscScalar *a;
84: INT_TYPE *ia, *ja;
86: /* Number of non-zero elements */
87: INT_TYPE nz;
89: /* Row permutaton vector*/
90: INT_TYPE *perm;
92: /* Define is matrix preserve sparce structure. */
93: MatStructure matstruc;
95: PetscErrorCode (*ConvertToTriples)(Mat, MatReuse, PetscInt*, PetscInt**, PetscInt**, PetscScalar**);
97: /* True if mkl_cpardiso function have been used. */
98: PetscBool CleanUp;
99: } Mat_MKL_CPARDISO;
101: /*
102: * Copy the elements of matrix A.
103: * Input:
104: * - Mat A: MATSEQAIJ matrix
105: * - int shift: matrix index.
106: * - 0 for c representation
107: * - 1 for fortran representation
108: * - MatReuse reuse:
109: * - MAT_INITIAL_MATRIX: Create a new aij representation
110: * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
111: * Output:
112: * - int *nnz: Number of nonzero-elements.
113: * - int **r pointer to i index
114: * - int **c pointer to j elements
115: * - MATRIXTYPE **v: Non-zero elements
116: */
117: PetscErrorCode MatCopy_seqaij_seqaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
118: {
119: Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data;
122: *v=aa->a;
123: if (reuse == MAT_INITIAL_MATRIX) {
124: *r = (INT_TYPE*)aa->i;
125: *c = (INT_TYPE*)aa->j;
126: *nnz = aa->nz;
127: }
128: return(0);
129: }
131: PetscErrorCode MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
132: {
133: const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
134: PetscErrorCode ierr;
135: PetscInt rstart,nz,i,j,countA,countB;
136: PetscInt *row,*col;
137: const PetscScalar *av, *bv;
138: PetscScalar *val;
139: Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data;
140: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(mat->A)->data;
141: Mat_SeqAIJ *bb = (Mat_SeqAIJ*)(mat->B)->data;
142: PetscInt colA_start,jB,jcol;
145: ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart;
146: av=aa->a; bv=bb->a;
148: garray = mat->garray;
150: if (reuse == MAT_INITIAL_MATRIX) {
151: nz = aa->nz + bb->nz;
152: *nnz = nz;
153: PetscMalloc3(m+1,&row,nz,&col,nz,&val);
154: *r = row; *c = col; *v = val;
155: } else {
156: row = *r; col = *c; val = *v;
157: }
159: nz = 0;
160: for (i=0; i<m; i++) {
161: row[i] = nz;
162: countA = ai[i+1] - ai[i];
163: countB = bi[i+1] - bi[i];
164: ajj = aj + ai[i]; /* ptr to the beginning of this row */
165: bjj = bj + bi[i];
167: /* B part, smaller col index */
168: colA_start = rstart + ajj[0]; /* the smallest global col index of A */
169: jB = 0;
170: for (j=0; j<countB; j++) {
171: jcol = garray[bjj[j]];
172: if (jcol > colA_start) break;
173: col[nz] = jcol;
174: val[nz++] = *bv++;
175: }
176: jB = j;
178: /* A part */
179: for (j=0; j<countA; j++) {
180: col[nz] = rstart + ajj[j];
181: val[nz++] = *av++;
182: }
184: /* B part, larger col index */
185: for (j=jB; j<countB; j++) {
186: col[nz] = garray[bjj[j]];
187: val[nz++] = *bv++;
188: }
189: }
190: row[m] = nz;
192: return(0);
193: }
195: PetscErrorCode MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
196: {
197: const PetscInt *ai, *aj, *bi, *bj,*garray,bs=A->rmap->bs,bs2=bs*bs,m=A->rmap->n/bs,*ajj,*bjj;
198: PetscErrorCode ierr;
199: PetscInt rstart,nz,i,j,countA,countB;
200: PetscInt *row,*col;
201: const PetscScalar *av, *bv;
202: PetscScalar *val;
203: Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data;
204: Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)(mat->A)->data;
205: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data;
206: PetscInt colA_start,jB,jcol;
209: ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart/bs;
210: av=aa->a; bv=bb->a;
212: garray = mat->garray;
214: if (reuse == MAT_INITIAL_MATRIX) {
215: nz = aa->nz + bb->nz;
216: *nnz = nz;
217: PetscMalloc3(m+1,&row,nz,&col,nz*bs2,&val);
218: *r = row; *c = col; *v = val;
219: } else {
220: row = *r; col = *c; val = *v;
221: }
223: nz = 0;
224: for (i=0; i<m; i++) {
225: row[i] = nz+1;
226: countA = ai[i+1] - ai[i];
227: countB = bi[i+1] - bi[i];
228: ajj = aj + ai[i]; /* ptr to the beginning of this row */
229: bjj = bj + bi[i];
231: /* B part, smaller col index */
232: colA_start = rstart + (countA > 0 ? ajj[0] : 0); /* the smallest global col index of A */
233: jB = 0;
234: for (j=0; j<countB; j++) {
235: jcol = garray[bjj[j]];
236: if (jcol > colA_start) break;
237: col[nz++] = jcol + 1;
238: }
239: jB = j;
240: PetscArraycpy(val,bv,jB*bs2);
241: val += jB*bs2;
242: bv += jB*bs2;
244: /* A part */
245: for (j=0; j<countA; j++) col[nz++] = rstart + ajj[j] + 1;
246: PetscArraycpy(val,av,countA*bs2);
247: val += countA*bs2;
248: av += countA*bs2;
250: /* B part, larger col index */
251: for (j=jB; j<countB; j++) col[nz++] = garray[bjj[j]] + 1;
252: PetscArraycpy(val,bv,(countB-jB)*bs2);
253: val += (countB-jB)*bs2;
254: bv += (countB-jB)*bs2;
255: }
256: row[m] = nz+1;
258: return(0);
259: }
261: PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
262: {
263: const PetscInt *ai, *aj, *bi, *bj,*garray,bs=A->rmap->bs,bs2=bs*bs,m=A->rmap->n/bs,*ajj,*bjj;
264: PetscErrorCode ierr;
265: PetscInt rstart,nz,i,j,countA,countB;
266: PetscInt *row,*col;
267: const PetscScalar *av, *bv;
268: PetscScalar *val;
269: Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data;
270: Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)(mat->A)->data;
271: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data;
274: ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart/bs;
275: av=aa->a; bv=bb->a;
277: garray = mat->garray;
279: if (reuse == MAT_INITIAL_MATRIX) {
280: nz = aa->nz + bb->nz;
281: *nnz = nz;
282: PetscMalloc3(m+1,&row,nz,&col,nz*bs2,&val);
283: *r = row; *c = col; *v = val;
284: } else {
285: row = *r; col = *c; val = *v;
286: }
288: nz = 0;
289: for (i=0; i<m; i++) {
290: row[i] = nz+1;
291: countA = ai[i+1] - ai[i];
292: countB = bi[i+1] - bi[i];
293: ajj = aj + ai[i]; /* ptr to the beginning of this row */
294: bjj = bj + bi[i];
296: /* A part */
297: for (j=0; j<countA; j++) col[nz++] = rstart + ajj[j] + 1;
298: PetscArraycpy(val,av,countA*bs2);
299: val += countA*bs2;
300: av += countA*bs2;
302: /* B part, larger col index */
303: for (j=0; j<countB; j++) col[nz++] = garray[bjj[j]] + 1;
304: PetscArraycpy(val,bv,countB*bs2);
305: val += countB*bs2;
306: bv += countB*bs2;
307: }
308: row[m] = nz+1;
310: return(0);
311: }
313: /*
314: * Free memory for Mat_MKL_CPARDISO structure and pointers to objects.
315: */
316: PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A)
317: {
318: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
319: PetscErrorCode ierr;
322: /* Terminate instance, deallocate memories */
323: if (mat_mkl_cpardiso->CleanUp) {
324: mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
326: cluster_sparse_solver (
327: mat_mkl_cpardiso->pt,
328: &mat_mkl_cpardiso->maxfct,
329: &mat_mkl_cpardiso->mnum,
330: &mat_mkl_cpardiso->mtype,
331: &mat_mkl_cpardiso->phase,
332: &mat_mkl_cpardiso->n,
333: NULL,
334: NULL,
335: NULL,
336: mat_mkl_cpardiso->perm,
337: &mat_mkl_cpardiso->nrhs,
338: mat_mkl_cpardiso->iparm,
339: &mat_mkl_cpardiso->msglvl,
340: NULL,
341: NULL,
342: &mat_mkl_cpardiso->comm_mkl_cpardiso,
343: (PetscInt*)&mat_mkl_cpardiso->err);
344: }
346: if (mat_mkl_cpardiso->ConvertToTriples != MatCopy_seqaij_seqaij_MKL_CPARDISO) {
347: PetscFree3(mat_mkl_cpardiso->ia,mat_mkl_cpardiso->ja,mat_mkl_cpardiso->a);
348: }
349: MPI_Comm_free(&(mat_mkl_cpardiso->comm_mkl_cpardiso));
350: PetscFree(A->data);
352: /* clear composed functions */
353: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
354: PetscObjectComposeFunction((PetscObject)A,"MatMkl_CPardisoSetCntl_C",NULL);
355: return(0);
356: }
358: /*
359: * Computes Ax = b
360: */
361: PetscErrorCode MatSolve_MKL_CPARDISO(Mat A,Vec b,Vec x)
362: {
363: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
364: PetscErrorCode ierr;
365: PetscScalar *xarray;
366: const PetscScalar *barray;
369: mat_mkl_cpardiso->nrhs = 1;
370: VecGetArray(x,&xarray);
371: VecGetArrayRead(b,&barray);
373: /* solve phase */
374: /*-------------*/
375: mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
376: cluster_sparse_solver (
377: mat_mkl_cpardiso->pt,
378: &mat_mkl_cpardiso->maxfct,
379: &mat_mkl_cpardiso->mnum,
380: &mat_mkl_cpardiso->mtype,
381: &mat_mkl_cpardiso->phase,
382: &mat_mkl_cpardiso->n,
383: mat_mkl_cpardiso->a,
384: mat_mkl_cpardiso->ia,
385: mat_mkl_cpardiso->ja,
386: mat_mkl_cpardiso->perm,
387: &mat_mkl_cpardiso->nrhs,
388: mat_mkl_cpardiso->iparm,
389: &mat_mkl_cpardiso->msglvl,
390: (void*)barray,
391: (void*)xarray,
392: &mat_mkl_cpardiso->comm_mkl_cpardiso,
393: (PetscInt*)&mat_mkl_cpardiso->err);
395: 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));
397: VecRestoreArray(x,&xarray);
398: VecRestoreArrayRead(b,&barray);
399: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
400: return(0);
401: }
403: PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A,Vec b,Vec x)
404: {
405: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
406: PetscErrorCode ierr;
409: #if defined(PETSC_USE_COMPLEX)
410: mat_mkl_cpardiso->iparm[12 - 1] = 1;
411: #else
412: mat_mkl_cpardiso->iparm[12 - 1] = 2;
413: #endif
414: MatSolve_MKL_CPARDISO(A,b,x);
415: mat_mkl_cpardiso->iparm[12 - 1] = 0;
416: return(0);
417: }
419: PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A,Mat B,Mat X)
420: {
421: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
422: PetscErrorCode ierr;
423: PetscScalar *xarray;
424: const PetscScalar *barray;
427: MatGetSize(B,NULL,(PetscInt*)&mat_mkl_cpardiso->nrhs);
429: if (mat_mkl_cpardiso->nrhs > 0) {
430: MatDenseGetArrayRead(B,&barray);
431: MatDenseGetArray(X,&xarray);
433: /* solve phase */
434: /*-------------*/
435: mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
436: cluster_sparse_solver (
437: mat_mkl_cpardiso->pt,
438: &mat_mkl_cpardiso->maxfct,
439: &mat_mkl_cpardiso->mnum,
440: &mat_mkl_cpardiso->mtype,
441: &mat_mkl_cpardiso->phase,
442: &mat_mkl_cpardiso->n,
443: mat_mkl_cpardiso->a,
444: mat_mkl_cpardiso->ia,
445: mat_mkl_cpardiso->ja,
446: mat_mkl_cpardiso->perm,
447: &mat_mkl_cpardiso->nrhs,
448: mat_mkl_cpardiso->iparm,
449: &mat_mkl_cpardiso->msglvl,
450: (void*)barray,
451: (void*)xarray,
452: &mat_mkl_cpardiso->comm_mkl_cpardiso,
453: (PetscInt*)&mat_mkl_cpardiso->err);
454: 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));
455: MatDenseRestoreArrayRead(B,&barray);
456: MatDenseRestoreArray(X,&xarray);
458: }
459: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
460: return(0);
462: }
464: /*
465: * LU Decomposition
466: */
467: PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F,Mat A,const MatFactorInfo *info)
468: {
469: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(F)->data;
470: PetscErrorCode ierr;
473: mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
474: (*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_REUSE_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);
476: mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION;
477: cluster_sparse_solver (
478: mat_mkl_cpardiso->pt,
479: &mat_mkl_cpardiso->maxfct,
480: &mat_mkl_cpardiso->mnum,
481: &mat_mkl_cpardiso->mtype,
482: &mat_mkl_cpardiso->phase,
483: &mat_mkl_cpardiso->n,
484: mat_mkl_cpardiso->a,
485: mat_mkl_cpardiso->ia,
486: mat_mkl_cpardiso->ja,
487: mat_mkl_cpardiso->perm,
488: &mat_mkl_cpardiso->nrhs,
489: mat_mkl_cpardiso->iparm,
490: &mat_mkl_cpardiso->msglvl,
491: NULL,
492: NULL,
493: &mat_mkl_cpardiso->comm_mkl_cpardiso,
494: &mat_mkl_cpardiso->err);
495: 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));
497: mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
498: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
499: return(0);
500: }
502: /* Sets mkl_cpardiso options from the options database */
503: PetscErrorCode PetscSetMKL_CPARDISOFromOptions(Mat F, Mat A)
504: {
505: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
506: PetscErrorCode ierr;
507: PetscInt icntl,threads;
508: PetscBool flg;
511: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_CPARDISO Options","Mat");
512: PetscOptionsInt("-mat_mkl_cpardiso_65","Number of threads to use","None",threads,&threads,&flg);
513: if (flg) mkl_set_num_threads((int)threads);
515: 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);
516: if (flg) mat_mkl_cpardiso->maxfct = icntl;
518: PetscOptionsInt("-mat_mkl_cpardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_cpardiso->mnum,&icntl,&flg);
519: if (flg) mat_mkl_cpardiso->mnum = icntl;
521: PetscOptionsInt("-mat_mkl_cpardiso_68","Message level information","None",mat_mkl_cpardiso->msglvl,&icntl,&flg);
522: if (flg) mat_mkl_cpardiso->msglvl = icntl;
524: PetscOptionsInt("-mat_mkl_cpardiso_69","Defines the matrix type","None",mat_mkl_cpardiso->mtype,&icntl,&flg);
525: if (flg) {
526: mat_mkl_cpardiso->mtype = icntl;
527: #if defined(PETSC_USE_REAL_SINGLE)
528: mat_mkl_cpardiso->iparm[27] = 1;
529: #else
530: mat_mkl_cpardiso->iparm[27] = 0;
531: #endif
532: mat_mkl_cpardiso->iparm[34] = 1;
533: }
534: PetscOptionsInt("-mat_mkl_cpardiso_1","Use default values","None",mat_mkl_cpardiso->iparm[0],&icntl,&flg);
536: if (flg && icntl != 0) {
537: PetscOptionsInt("-mat_mkl_cpardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_cpardiso->iparm[1],&icntl,&flg);
538: if (flg) mat_mkl_cpardiso->iparm[1] = icntl;
540: PetscOptionsInt("-mat_mkl_cpardiso_4","Preconditioned CGS/CG","None",mat_mkl_cpardiso->iparm[3],&icntl,&flg);
541: if (flg) mat_mkl_cpardiso->iparm[3] = icntl;
543: PetscOptionsInt("-mat_mkl_cpardiso_5","User permutation","None",mat_mkl_cpardiso->iparm[4],&icntl,&flg);
544: if (flg) mat_mkl_cpardiso->iparm[4] = icntl;
546: PetscOptionsInt("-mat_mkl_cpardiso_6","Write solution on x","None",mat_mkl_cpardiso->iparm[5],&icntl,&flg);
547: if (flg) mat_mkl_cpardiso->iparm[5] = icntl;
549: PetscOptionsInt("-mat_mkl_cpardiso_8","Iterative refinement step","None",mat_mkl_cpardiso->iparm[7],&icntl,&flg);
550: if (flg) mat_mkl_cpardiso->iparm[7] = icntl;
552: PetscOptionsInt("-mat_mkl_cpardiso_10","Pivoting perturbation","None",mat_mkl_cpardiso->iparm[9],&icntl,&flg);
553: if (flg) mat_mkl_cpardiso->iparm[9] = icntl;
555: PetscOptionsInt("-mat_mkl_cpardiso_11","Scaling vectors","None",mat_mkl_cpardiso->iparm[10],&icntl,&flg);
556: if (flg) mat_mkl_cpardiso->iparm[10] = icntl;
558: PetscOptionsInt("-mat_mkl_cpardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_cpardiso->iparm[11],&icntl,&flg);
559: if (flg) mat_mkl_cpardiso->iparm[11] = icntl;
561: PetscOptionsInt("-mat_mkl_cpardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_cpardiso->iparm[12],&icntl,
562: &flg);
563: if (flg) mat_mkl_cpardiso->iparm[12] = icntl;
565: PetscOptionsInt("-mat_mkl_cpardiso_18","Numbers of non-zero elements","None",mat_mkl_cpardiso->iparm[17],&icntl,
566: &flg);
567: if (flg) mat_mkl_cpardiso->iparm[17] = icntl;
569: PetscOptionsInt("-mat_mkl_cpardiso_19","Report number of floating point operations","None",mat_mkl_cpardiso->iparm[18],&icntl,&flg);
570: if (flg) mat_mkl_cpardiso->iparm[18] = icntl;
572: PetscOptionsInt("-mat_mkl_cpardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_cpardiso->iparm[20],&icntl,&flg);
573: if (flg) mat_mkl_cpardiso->iparm[20] = icntl;
575: PetscOptionsInt("-mat_mkl_cpardiso_24","Parallel factorization control","None",mat_mkl_cpardiso->iparm[23],&icntl,&flg);
576: if (flg) mat_mkl_cpardiso->iparm[23] = icntl;
578: PetscOptionsInt("-mat_mkl_cpardiso_25","Parallel forward/backward solve control","None",mat_mkl_cpardiso->iparm[24],&icntl,&flg);
579: if (flg) mat_mkl_cpardiso->iparm[24] = icntl;
581: PetscOptionsInt("-mat_mkl_cpardiso_27","Matrix checker","None",mat_mkl_cpardiso->iparm[26],&icntl,&flg);
582: if (flg) mat_mkl_cpardiso->iparm[26] = icntl;
584: PetscOptionsInt("-mat_mkl_cpardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_cpardiso->iparm[30],&icntl,&flg);
585: if (flg) mat_mkl_cpardiso->iparm[30] = icntl;
587: PetscOptionsInt("-mat_mkl_cpardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_cpardiso->iparm[33],&icntl,&flg);
588: if (flg) mat_mkl_cpardiso->iparm[33] = icntl;
590: PetscOptionsInt("-mat_mkl_cpardiso_60","Intel MKL_CPARDISO mode","None",mat_mkl_cpardiso->iparm[59],&icntl,&flg);
591: if (flg) mat_mkl_cpardiso->iparm[59] = icntl;
592: }
594: PetscOptionsEnd();
595: return(0);
596: }
598: PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso)
599: {
600: PetscErrorCode ierr;
601: PetscInt bs;
602: PetscBool match;
603: PetscMPIInt size;
607: MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mat_mkl_cpardiso->comm_mkl_cpardiso));
608: MPI_Comm_size(mat_mkl_cpardiso->comm_mkl_cpardiso, &size);
610: mat_mkl_cpardiso->CleanUp = PETSC_FALSE;
611: mat_mkl_cpardiso->maxfct = 1;
612: mat_mkl_cpardiso->mnum = 1;
613: mat_mkl_cpardiso->n = A->rmap->N;
614: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
615: mat_mkl_cpardiso->msglvl = 0;
616: mat_mkl_cpardiso->nrhs = 1;
617: mat_mkl_cpardiso->err = 0;
618: mat_mkl_cpardiso->phase = -1;
619: #if defined(PETSC_USE_COMPLEX)
620: mat_mkl_cpardiso->mtype = 13;
621: #else
622: mat_mkl_cpardiso->mtype = 11;
623: #endif
625: #if defined(PETSC_USE_REAL_SINGLE)
626: mat_mkl_cpardiso->iparm[27] = 1;
627: #else
628: mat_mkl_cpardiso->iparm[27] = 0;
629: #endif
631: mat_mkl_cpardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */
632: mat_mkl_cpardiso->iparm[ 1] = 2; /* Use METIS for fill-in reordering */
633: mat_mkl_cpardiso->iparm[ 5] = 0; /* Write solution into x */
634: mat_mkl_cpardiso->iparm[ 7] = 2; /* Max number of iterative refinement steps */
635: mat_mkl_cpardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
636: mat_mkl_cpardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
637: mat_mkl_cpardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
638: mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
639: mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
640: mat_mkl_cpardiso->iparm[26] = 1; /* Check input data for correctness */
642: mat_mkl_cpardiso->iparm[39] = 0;
643: if (size > 1) {
644: mat_mkl_cpardiso->iparm[39] = 2;
645: mat_mkl_cpardiso->iparm[40] = A->rmap->rstart;
646: mat_mkl_cpardiso->iparm[41] = A->rmap->rend-1;
647: }
648: PetscObjectTypeCompareAny((PetscObject)A,&match,MATMPIBAIJ,MATMPISBAIJ,"");
649: if (match) {
650: MatGetBlockSize(A,&bs);
651: mat_mkl_cpardiso->iparm[36] = bs;
652: mat_mkl_cpardiso->iparm[40] /= bs;
653: mat_mkl_cpardiso->iparm[41] /= bs;
654: mat_mkl_cpardiso->iparm[40]++;
655: mat_mkl_cpardiso->iparm[41]++;
656: mat_mkl_cpardiso->iparm[34] = 0; /* Fortran style */
657: } else {
658: mat_mkl_cpardiso->iparm[34] = 1; /* C style */
659: }
661: mat_mkl_cpardiso->perm = 0;
662: return(0);
663: }
665: /*
666: * Symbolic decomposition. Mkl_Pardiso analysis phase.
667: */
668: PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
669: {
670: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
671: PetscErrorCode ierr;
674: mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
676: /* Set MKL_CPARDISO options from the options database */
677: PetscSetMKL_CPARDISOFromOptions(F,A);
678: (*mat_mkl_cpardiso->ConvertToTriples)(A,MAT_INITIAL_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);
680: mat_mkl_cpardiso->n = A->rmap->N;
681: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
683: /* analysis phase */
684: /*----------------*/
685: mat_mkl_cpardiso->phase = JOB_ANALYSIS;
687: cluster_sparse_solver (
688: mat_mkl_cpardiso->pt,
689: &mat_mkl_cpardiso->maxfct,
690: &mat_mkl_cpardiso->mnum,
691: &mat_mkl_cpardiso->mtype,
692: &mat_mkl_cpardiso->phase,
693: &mat_mkl_cpardiso->n,
694: mat_mkl_cpardiso->a,
695: mat_mkl_cpardiso->ia,
696: mat_mkl_cpardiso->ja,
697: mat_mkl_cpardiso->perm,
698: &mat_mkl_cpardiso->nrhs,
699: mat_mkl_cpardiso->iparm,
700: &mat_mkl_cpardiso->msglvl,
701: NULL,
702: NULL,
703: &mat_mkl_cpardiso->comm_mkl_cpardiso,
704: (PetscInt*)&mat_mkl_cpardiso->err);
706: 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));
708: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
709: F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO;
710: F->ops->solve = MatSolve_MKL_CPARDISO;
711: F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO;
712: F->ops->matsolve = MatMatSolve_MKL_CPARDISO;
713: return(0);
714: }
716: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS perm,const MatFactorInfo *info)
717: {
718: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
719: PetscErrorCode ierr;
722: mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
725: /* Set MKL_CPARDISO options from the options database */
726: PetscSetMKL_CPARDISOFromOptions(F,A);
727: (*mat_mkl_cpardiso->ConvertToTriples)(A,MAT_INITIAL_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);
729: mat_mkl_cpardiso->n = A->rmap->N;
730: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
731: #if defined(PETSC_USE_COMPLEX)
732: SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead",((PetscObject)A)->type_name);
733: #endif
734: if (A->spd_set && A->spd) mat_mkl_cpardiso->mtype = 2;
735: else mat_mkl_cpardiso->mtype = -2;
737: /* analysis phase */
738: /*----------------*/
739: mat_mkl_cpardiso->phase = JOB_ANALYSIS;
741: cluster_sparse_solver (
742: mat_mkl_cpardiso->pt,
743: &mat_mkl_cpardiso->maxfct,
744: &mat_mkl_cpardiso->mnum,
745: &mat_mkl_cpardiso->mtype,
746: &mat_mkl_cpardiso->phase,
747: &mat_mkl_cpardiso->n,
748: mat_mkl_cpardiso->a,
749: mat_mkl_cpardiso->ia,
750: mat_mkl_cpardiso->ja,
751: mat_mkl_cpardiso->perm,
752: &mat_mkl_cpardiso->nrhs,
753: mat_mkl_cpardiso->iparm,
754: &mat_mkl_cpardiso->msglvl,
755: NULL,
756: NULL,
757: &mat_mkl_cpardiso->comm_mkl_cpardiso,
758: (PetscInt*)&mat_mkl_cpardiso->err);
760: 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));
762: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
763: F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_CPARDISO;
764: F->ops->solve = MatSolve_MKL_CPARDISO;
765: F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO;
766: F->ops->matsolve = MatMatSolve_MKL_CPARDISO;
767: return(0);
768: }
770: PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer)
771: {
772: PetscErrorCode ierr;
773: PetscBool iascii;
774: PetscViewerFormat format;
775: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
776: PetscInt i;
779: /* check if matrix is mkl_cpardiso type */
780: if (A->ops->solve != MatSolve_MKL_CPARDISO) return(0);
782: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
783: if (iascii) {
784: PetscViewerGetFormat(viewer,&format);
785: if (format == PETSC_VIEWER_ASCII_INFO) {
786: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO run parameters:\n");
787: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO phase: %d \n",mat_mkl_cpardiso->phase);
788: for (i = 1; i <= 64; i++) {
789: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO iparm[%d]: %d \n",i, mat_mkl_cpardiso->iparm[i - 1]);
790: }
791: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO maxfct: %d \n", mat_mkl_cpardiso->maxfct);
792: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mnum: %d \n", mat_mkl_cpardiso->mnum);
793: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mtype: %d \n", mat_mkl_cpardiso->mtype);
794: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO n: %d \n", mat_mkl_cpardiso->n);
795: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO nrhs: %d \n", mat_mkl_cpardiso->nrhs);
796: PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO msglvl: %d \n", mat_mkl_cpardiso->msglvl);
797: }
798: }
799: return(0);
800: }
802: PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info)
803: {
804: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
807: info->block_size = 1.0;
808: info->nz_allocated = mat_mkl_cpardiso->nz + 0.0;
809: info->nz_unneeded = 0.0;
810: info->assemblies = 0.0;
811: info->mallocs = 0.0;
812: info->memory = 0.0;
813: info->fill_ratio_given = 0;
814: info->fill_ratio_needed = 0;
815: info->factor_mallocs = 0;
816: return(0);
817: }
819: PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F,PetscInt icntl,PetscInt ival)
820: {
821: Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)F->data;
824: if (icntl <= 64) {
825: mat_mkl_cpardiso->iparm[icntl - 1] = ival;
826: } else {
827: if (icntl == 65) mkl_set_num_threads((int)ival);
828: else if (icntl == 66) mat_mkl_cpardiso->maxfct = ival;
829: else if (icntl == 67) mat_mkl_cpardiso->mnum = ival;
830: else if (icntl == 68) mat_mkl_cpardiso->msglvl = ival;
831: else if (icntl == 69) {
832: mat_mkl_cpardiso->mtype = ival;
833: #if defined(PETSC_USE_REAL_SINGLE)
834: mat_mkl_cpardiso->iparm[27] = 1;
835: #else
836: mat_mkl_cpardiso->iparm[27] = 0;
837: #endif
838: mat_mkl_cpardiso->iparm[34] = 1;
839: }
840: }
841: return(0);
842: }
844: /*@
845: MatMkl_CPardisoSetCntl - Set Mkl_Pardiso parameters
847: Logically Collective on Mat
849: Input Parameters:
850: + F - the factored matrix obtained by calling MatGetFactor()
851: . icntl - index of Mkl_Pardiso parameter
852: - ival - value of Mkl_Pardiso parameter
854: Options Database:
855: . -mat_mkl_cpardiso_<icntl> <ival>
857: Level: Intermediate
859: Notes:
860: 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
861: database approach when working with TS, SNES, or KSP.
863: References:
864: . Mkl_Pardiso Users' Guide
866: .seealso: MatGetFactor()
867: @*/
868: PetscErrorCode MatMkl_CPardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
869: {
873: PetscTryMethod(F,"MatMkl_CPardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
874: return(0);
875: }
877: static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type)
878: {
880: *type = MATSOLVERMKL_CPARDISO;
881: return(0);
882: }
884: /* MatGetFactor for MPI AIJ matrices */
885: static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A,MatFactorType ftype,Mat *F)
886: {
887: Mat B;
888: PetscErrorCode ierr;
889: Mat_MKL_CPARDISO *mat_mkl_cpardiso;
890: PetscBool isSeqAIJ,isMPIBAIJ,isMPISBAIJ;
893: /* Create the factorization matrix */
895: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
896: PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&isMPIBAIJ);
897: PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMPISBAIJ);
898: MatCreate(PetscObjectComm((PetscObject)A),&B);
899: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
900: PetscStrallocpy("mkl_cpardiso",&((PetscObject)B)->type_name);
901: MatSetUp(B);
903: PetscNewLog(B,&mat_mkl_cpardiso);
905: if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO;
906: else if (isMPIBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO;
907: else if (isMPISBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO;
908: else mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO;
910: if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO;
911: else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_CPARDISO;
912: B->ops->destroy = MatDestroy_MKL_CPARDISO;
914: B->ops->view = MatView_MKL_CPARDISO;
915: B->ops->getinfo = MatGetInfo_MKL_CPARDISO;
917: B->factortype = ftype;
918: B->assembled = PETSC_TRUE; /* required by -ksp_view */
920: B->data = mat_mkl_cpardiso;
922: /* set solvertype */
923: PetscFree(B->solvertype);
924: PetscStrallocpy(MATSOLVERMKL_CPARDISO,&B->solvertype);
926: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_cpardiso);
927: PetscObjectComposeFunction((PetscObject)B,"MatMkl_CPardisoSetCntl_C",MatMkl_CPardisoSetCntl_MKL_CPARDISO);
928: PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso);
930: *F = B;
931: return(0);
932: }
934: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void)
935: {
937:
939: MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
940: MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
941: MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIBAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
942: MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPISBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_mpiaij_mkl_cpardiso);
943: return(0);
944: }