Actual source code: mkl_pardiso.c
1: #include <../src/mat/impls/aij/seq/aij.h>
2: #include <../src/mat/impls/sbaij/seq/sbaij.h>
3: #include <../src/mat/impls/dense/seq/dense.h>
5: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
6: #define MKL_ILP64
7: #endif
8: #include <mkl_pardiso.h>
10: PETSC_EXTERN void PetscSetMKL_PARDISOThreads(int);
12: /*
13: * Possible mkl_pardiso phases that controls the execution of the solver.
14: * For more information check mkl_pardiso 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
30: #if defined(PETSC_USE_64BIT_INDICES)
31: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
32: #define INT_TYPE long long int
33: #define MKL_PARDISO pardiso
34: #define MKL_PARDISO_INIT pardisoinit
35: #else
36: /* this is the case where the MKL BLAS/LAPACK are 32 bit integers but the 64 bit integer version of
37: of Pardiso code is used; hence the need for the 64 below*/
38: #define INT_TYPE long long int
39: #define MKL_PARDISO pardiso_64
40: #define MKL_PARDISO_INIT pardiso_64init
41: void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm [])
42: {
43: int iparm_copy[IPARM_SIZE], mtype_copy, i;
45: mtype_copy = *mtype;
46: pardisoinit(pt, &mtype_copy, iparm_copy);
47: for (i=0; i<IPARM_SIZE; i++) iparm[i] = iparm_copy[i];
48: }
49: #endif
50: #else
51: #define INT_TYPE int
52: #define MKL_PARDISO pardiso
53: #define MKL_PARDISO_INIT pardisoinit
54: #endif
56: /*
57: * Internal data structure.
58: * For more information check mkl_pardiso manual.
59: */
60: typedef struct {
62: /* Configuration vector*/
63: INT_TYPE iparm[IPARM_SIZE];
65: /*
66: * Internal mkl_pardiso memory location.
67: * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
68: */
69: void *pt[IPARM_SIZE];
71: /* Basic mkl_pardiso info*/
72: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
74: /* Matrix structure*/
75: void *a;
76: INT_TYPE *ia, *ja;
78: /* Number of non-zero elements*/
79: INT_TYPE nz;
81: /* Row permutaton vector*/
82: INT_TYPE *perm;
84: /* Define if matrix preserves sparse structure.*/
85: MatStructure matstruc;
87: PetscBool needsym;
88: PetscBool freeaij;
90: /* Schur complement */
91: PetscScalar *schur;
92: PetscInt schur_size;
93: PetscInt *schur_idxs;
94: PetscScalar *schur_work;
95: PetscBLASInt schur_work_size;
96: PetscBool solve_interior;
98: /* True if mkl_pardiso function have been used.*/
99: PetscBool CleanUp;
101: /* Conversion to a format suitable for MKL */
102: PetscErrorCode (*Convert)(Mat, PetscBool, MatReuse, PetscBool*, INT_TYPE*, INT_TYPE**, INT_TYPE**, PetscScalar**);
103: } Mat_MKL_PARDISO;
105: PetscErrorCode MatMKLPardiso_Convert_seqsbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
106: {
107: Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data;
108: PetscInt bs = A->rmap->bs,i;
112: if (!sym) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
113: *v = aa->a;
114: if (bs == 1) { /* already in the correct format */
115: /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
116: *r = (INT_TYPE*)aa->i;
117: *c = (INT_TYPE*)aa->j;
118: *nnz = (INT_TYPE)aa->nz;
119: *free = PETSC_FALSE;
120: } else if (reuse == MAT_INITIAL_MATRIX) {
121: PetscInt m = A->rmap->n,nz = aa->nz;
122: PetscInt *row,*col;
123: PetscMalloc2(m+1,&row,nz,&col);
124: for (i=0; i<m+1; i++) {
125: row[i] = aa->i[i]+1;
126: }
127: for (i=0; i<nz; i++) {
128: col[i] = aa->j[i]+1;
129: }
130: *r = (INT_TYPE*)row;
131: *c = (INT_TYPE*)col;
132: *nnz = (INT_TYPE)nz;
133: *free = PETSC_TRUE;
134: }
135: return(0);
136: }
138: PetscErrorCode MatMKLPardiso_Convert_seqbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
139: {
140: Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)A->data;
141: PetscInt bs = A->rmap->bs,i;
145: if (!sym) {
146: *v = aa->a;
147: if (bs == 1) { /* already in the correct format */
148: /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
149: *r = (INT_TYPE*)aa->i;
150: *c = (INT_TYPE*)aa->j;
151: *nnz = (INT_TYPE)aa->nz;
152: *free = PETSC_FALSE;
153: return(0);
154: } else if (reuse == MAT_INITIAL_MATRIX) {
155: PetscInt m = A->rmap->n,nz = aa->nz;
156: PetscInt *row,*col;
157: PetscMalloc2(m+1,&row,nz,&col);
158: for (i=0; i<m+1; i++) {
159: row[i] = aa->i[i]+1;
160: }
161: for (i=0; i<nz; i++) {
162: col[i] = aa->j[i]+1;
163: }
164: *r = (INT_TYPE*)row;
165: *c = (INT_TYPE*)col;
166: *nnz = (INT_TYPE)nz;
167: }
168: *free = PETSC_TRUE;
169: } else {
170: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
171: }
172: return(0);
173: }
175: PetscErrorCode MatMKLPardiso_Convert_seqaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
176: {
177: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data;
178: PetscScalar *aav;
182: MatSeqAIJGetArrayRead(A,(const PetscScalar**)&aav);
183: if (!sym) { /* already in the correct format */
184: *v = aav;
185: *r = (INT_TYPE*)aa->i;
186: *c = (INT_TYPE*)aa->j;
187: *nnz = (INT_TYPE)aa->nz;
188: *free = PETSC_FALSE;
189: } else if (reuse == MAT_INITIAL_MATRIX) { /* need to get the triangular part */
190: PetscScalar *vals,*vv;
191: PetscInt *row,*col,*jj;
192: PetscInt m = A->rmap->n,nz,i;
194: nz = 0;
195: for (i=0; i<m; i++) nz += aa->i[i+1] - aa->diag[i];
196: PetscMalloc2(m+1,&row,nz,&col);
197: PetscMalloc1(nz,&vals);
198: jj = col;
199: vv = vals;
201: row[0] = 0;
202: for (i=0; i<m; i++) {
203: PetscInt *aj = aa->j + aa->diag[i];
204: PetscScalar *av = aav + aa->diag[i];
205: PetscInt rl = aa->i[i+1] - aa->diag[i],j;
207: for (j=0; j<rl; j++) {
208: *jj = *aj; jj++; aj++;
209: *vv = *av; vv++; av++;
210: }
211: row[i+1] = row[i] + rl;
212: }
213: *v = vals;
214: *r = (INT_TYPE*)row;
215: *c = (INT_TYPE*)col;
216: *nnz = (INT_TYPE)nz;
217: *free = PETSC_TRUE;
218: } else {
219: PetscScalar *vv;
220: PetscInt m = A->rmap->n,i;
222: vv = *v;
223: for (i=0; i<m; i++) {
224: PetscScalar *av = aav + aa->diag[i];
225: PetscInt rl = aa->i[i+1] - aa->diag[i],j;
226: for (j=0; j<rl; j++) {
227: *vv = *av; vv++; av++;
228: }
229: }
230: *free = PETSC_TRUE;
231: }
232: MatSeqAIJRestoreArrayRead(A,(const PetscScalar**)&aav);
233: return(0);
234: }
236: static PetscErrorCode MatMKLPardisoSolveSchur_Private(Mat F, PetscScalar *B, PetscScalar *X)
237: {
238: Mat_MKL_PARDISO *mpardiso = (Mat_MKL_PARDISO*)F->data;
239: Mat S,Xmat,Bmat;
240: MatFactorSchurStatus schurstatus;
241: PetscErrorCode ierr;
244: MatFactorGetSchurComplement(F,&S,&schurstatus);
245: if (X == B && schurstatus == MAT_FACTOR_SCHUR_INVERTED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"X and B cannot point to the same address");
246: MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,B,&Bmat);
247: MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,X,&Xmat);
248: MatSetType(Bmat,((PetscObject)S)->type_name);
249: MatSetType(Xmat,((PetscObject)S)->type_name);
250: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
251: MatBindToCPU(Xmat,S->boundtocpu);
252: MatBindToCPU(Bmat,S->boundtocpu);
253: #endif
255: #if defined(PETSC_USE_COMPLEX)
256: if (mpardiso->iparm[12-1] == 1) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Hermitian solve not implemented yet");
257: #endif
259: switch (schurstatus) {
260: case MAT_FACTOR_SCHUR_FACTORED:
261: if (!mpardiso->iparm[12-1]) {
262: MatMatSolve(S,Bmat,Xmat);
263: } else { /* transpose solve */
264: MatMatSolveTranspose(S,Bmat,Xmat);
265: }
266: break;
267: case MAT_FACTOR_SCHUR_INVERTED:
268: MatProductCreateWithMat(S,Bmat,NULL,Xmat);
269: if (!mpardiso->iparm[12-1]) {
270: MatProductSetType(Xmat,MATPRODUCT_AB);
271: } else { /* transpose solve */
272: MatProductSetType(Xmat,MATPRODUCT_AtB);
273: }
274: MatProductSetFromOptions(Xmat);
275: MatProductSymbolic(Xmat);
276: MatProductNumeric(Xmat);
277: MatProductClear(Xmat);
278: break;
279: default:
280: SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status);
281: break;
282: }
283: MatFactorRestoreSchurComplement(F,&S,schurstatus);
284: MatDestroy(&Bmat);
285: MatDestroy(&Xmat);
286: return(0);
287: }
289: PetscErrorCode MatFactorSetSchurIS_MKL_PARDISO(Mat F, IS is)
290: {
291: Mat_MKL_PARDISO *mpardiso = (Mat_MKL_PARDISO*)F->data;
292: const PetscScalar *arr;
293: const PetscInt *idxs;
294: PetscInt size,i;
295: PetscMPIInt csize;
296: PetscBool sorted;
297: PetscErrorCode ierr;
300: MPI_Comm_size(PetscObjectComm((PetscObject)F),&csize);
301: if (csize > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MKL_PARDISO parallel Schur complements not yet supported from PETSc");
302: ISSorted(is,&sorted);
303: if (!sorted) {
304: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS for MKL_PARDISO Schur complements needs to be sorted");
305: }
306: ISGetLocalSize(is,&size);
307: PetscFree(mpardiso->schur_work);
308: PetscBLASIntCast(PetscMax(mpardiso->n,2*size),&mpardiso->schur_work_size);
309: PetscMalloc1(mpardiso->schur_work_size,&mpardiso->schur_work);
310: MatDestroy(&F->schur);
311: MatCreateSeqDense(PETSC_COMM_SELF,size,size,NULL,&F->schur);
312: MatDenseGetArrayRead(F->schur,&arr);
313: mpardiso->schur = (PetscScalar*)arr;
314: mpardiso->schur_size = size;
315: MatDenseRestoreArrayRead(F->schur,&arr);
316: if (mpardiso->mtype == 2) {
317: MatSetOption(F->schur,MAT_SPD,PETSC_TRUE);
318: }
320: PetscFree(mpardiso->schur_idxs);
321: PetscMalloc1(size,&mpardiso->schur_idxs);
322: PetscArrayzero(mpardiso->perm,mpardiso->n);
323: ISGetIndices(is,&idxs);
324: PetscArraycpy(mpardiso->schur_idxs,idxs,size);
325: for (i=0;i<size;i++) mpardiso->perm[idxs[i]] = 1;
326: ISRestoreIndices(is,&idxs);
327: if (size) { /* turn on Schur switch if the set of indices is not empty */
328: mpardiso->iparm[36-1] = 2;
329: }
330: return(0);
331: }
333: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
334: {
335: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
336: PetscErrorCode ierr;
339: if (mat_mkl_pardiso->CleanUp) {
340: mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
342: MKL_PARDISO (mat_mkl_pardiso->pt,
343: &mat_mkl_pardiso->maxfct,
344: &mat_mkl_pardiso->mnum,
345: &mat_mkl_pardiso->mtype,
346: &mat_mkl_pardiso->phase,
347: &mat_mkl_pardiso->n,
348: NULL,
349: NULL,
350: NULL,
351: NULL,
352: &mat_mkl_pardiso->nrhs,
353: mat_mkl_pardiso->iparm,
354: &mat_mkl_pardiso->msglvl,
355: NULL,
356: NULL,
357: &mat_mkl_pardiso->err);
358: }
359: PetscFree(mat_mkl_pardiso->perm);
360: PetscFree(mat_mkl_pardiso->schur_work);
361: PetscFree(mat_mkl_pardiso->schur_idxs);
362: if (mat_mkl_pardiso->freeaij) {
363: PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
364: if (mat_mkl_pardiso->iparm[34] == 1) {
365: PetscFree(mat_mkl_pardiso->a);
366: }
367: }
368: PetscFree(A->data);
370: /* clear composed functions */
371: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
372: PetscObjectComposeFunction((PetscObject)A,"MatFactorSetSchurIS_C",NULL);
373: PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);
374: return(0);
375: }
377: static PetscErrorCode MatMKLPardisoScatterSchur_Private(Mat_MKL_PARDISO *mpardiso, PetscScalar *whole, PetscScalar *schur, PetscBool reduce)
378: {
380: if (reduce) { /* data given for the whole matrix */
381: PetscInt i,m=0,p=0;
382: for (i=0;i<mpardiso->nrhs;i++) {
383: PetscInt j;
384: for (j=0;j<mpardiso->schur_size;j++) {
385: schur[p+j] = whole[m+mpardiso->schur_idxs[j]];
386: }
387: m += mpardiso->n;
388: p += mpardiso->schur_size;
389: }
390: } else { /* from Schur to whole */
391: PetscInt i,m=0,p=0;
392: for (i=0;i<mpardiso->nrhs;i++) {
393: PetscInt j;
394: for (j=0;j<mpardiso->schur_size;j++) {
395: whole[m+mpardiso->schur_idxs[j]] = schur[p+j];
396: }
397: m += mpardiso->n;
398: p += mpardiso->schur_size;
399: }
400: }
401: return(0);
402: }
404: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x)
405: {
406: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
407: PetscErrorCode ierr;
408: PetscScalar *xarray;
409: const PetscScalar *barray;
412: mat_mkl_pardiso->nrhs = 1;
413: VecGetArrayWrite(x,&xarray);
414: VecGetArrayRead(b,&barray);
416: if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
417: else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
419: if (barray == xarray) { /* if the two vectors share the same memory */
420: PetscScalar *work;
421: if (!mat_mkl_pardiso->schur_work) {
422: PetscMalloc1(mat_mkl_pardiso->n,&work);
423: } else {
424: work = mat_mkl_pardiso->schur_work;
425: }
426: mat_mkl_pardiso->iparm[6-1] = 1;
427: MKL_PARDISO (mat_mkl_pardiso->pt,
428: &mat_mkl_pardiso->maxfct,
429: &mat_mkl_pardiso->mnum,
430: &mat_mkl_pardiso->mtype,
431: &mat_mkl_pardiso->phase,
432: &mat_mkl_pardiso->n,
433: mat_mkl_pardiso->a,
434: mat_mkl_pardiso->ia,
435: mat_mkl_pardiso->ja,
436: NULL,
437: &mat_mkl_pardiso->nrhs,
438: mat_mkl_pardiso->iparm,
439: &mat_mkl_pardiso->msglvl,
440: (void*)xarray,
441: (void*)work,
442: &mat_mkl_pardiso->err);
443: if (!mat_mkl_pardiso->schur_work) {
444: PetscFree(work);
445: }
446: } else {
447: mat_mkl_pardiso->iparm[6-1] = 0;
448: MKL_PARDISO (mat_mkl_pardiso->pt,
449: &mat_mkl_pardiso->maxfct,
450: &mat_mkl_pardiso->mnum,
451: &mat_mkl_pardiso->mtype,
452: &mat_mkl_pardiso->phase,
453: &mat_mkl_pardiso->n,
454: mat_mkl_pardiso->a,
455: mat_mkl_pardiso->ia,
456: mat_mkl_pardiso->ja,
457: mat_mkl_pardiso->perm,
458: &mat_mkl_pardiso->nrhs,
459: mat_mkl_pardiso->iparm,
460: &mat_mkl_pardiso->msglvl,
461: (void*)barray,
462: (void*)xarray,
463: &mat_mkl_pardiso->err);
464: }
465: VecRestoreArrayRead(b,&barray);
467: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
469: if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
470: if (!mat_mkl_pardiso->solve_interior) {
471: PetscInt shift = mat_mkl_pardiso->schur_size;
473: MatFactorFactorizeSchurComplement(A);
474: /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
475: if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;
477: /* solve Schur complement */
478: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
479: MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
480: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
481: } else { /* if we are solving for the interior problem, any value in barray[schur] forward-substituted to xarray[schur] will be neglected */
482: PetscInt i;
483: for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
484: xarray[mat_mkl_pardiso->schur_idxs[i]] = 0.;
485: }
486: }
488: /* expansion phase */
489: mat_mkl_pardiso->iparm[6-1] = 1;
490: mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
491: MKL_PARDISO (mat_mkl_pardiso->pt,
492: &mat_mkl_pardiso->maxfct,
493: &mat_mkl_pardiso->mnum,
494: &mat_mkl_pardiso->mtype,
495: &mat_mkl_pardiso->phase,
496: &mat_mkl_pardiso->n,
497: mat_mkl_pardiso->a,
498: mat_mkl_pardiso->ia,
499: mat_mkl_pardiso->ja,
500: mat_mkl_pardiso->perm,
501: &mat_mkl_pardiso->nrhs,
502: mat_mkl_pardiso->iparm,
503: &mat_mkl_pardiso->msglvl,
504: (void*)xarray,
505: (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
506: &mat_mkl_pardiso->err);
508: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
509: mat_mkl_pardiso->iparm[6-1] = 0;
510: }
511: VecRestoreArrayWrite(x,&xarray);
512: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
513: return(0);
514: }
516: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x)
517: {
518: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
519: PetscInt oiparm12;
520: PetscErrorCode ierr;
523: oiparm12 = mat_mkl_pardiso->iparm[12 - 1];
524: mat_mkl_pardiso->iparm[12 - 1] = 2;
525: MatSolve_MKL_PARDISO(A,b,x);
526: mat_mkl_pardiso->iparm[12 - 1] = oiparm12;
527: return(0);
528: }
530: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X)
531: {
532: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->data;
533: PetscErrorCode ierr;
534: const PetscScalar *barray;
535: PetscScalar *xarray;
536: PetscBool flg;
539: PetscObjectBaseTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
540: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
541: if (X != B) {
542: PetscObjectBaseTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
543: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");
544: }
546: MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);
548: if (mat_mkl_pardiso->nrhs > 0) {
549: MatDenseGetArrayRead(B,&barray);
550: MatDenseGetArrayWrite(X,&xarray);
552: if (barray == xarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"B and X cannot share the same memory location");
553: if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
554: else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
556: MKL_PARDISO (mat_mkl_pardiso->pt,
557: &mat_mkl_pardiso->maxfct,
558: &mat_mkl_pardiso->mnum,
559: &mat_mkl_pardiso->mtype,
560: &mat_mkl_pardiso->phase,
561: &mat_mkl_pardiso->n,
562: mat_mkl_pardiso->a,
563: mat_mkl_pardiso->ia,
564: mat_mkl_pardiso->ja,
565: mat_mkl_pardiso->perm,
566: &mat_mkl_pardiso->nrhs,
567: mat_mkl_pardiso->iparm,
568: &mat_mkl_pardiso->msglvl,
569: (void*)barray,
570: (void*)xarray,
571: &mat_mkl_pardiso->err);
572: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
574: MatDenseRestoreArrayRead(B,&barray);
575: if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
576: PetscScalar *o_schur_work = NULL;
578: /* solve Schur complement */
579: if (!mat_mkl_pardiso->solve_interior) {
580: PetscInt shift = mat_mkl_pardiso->schur_size*mat_mkl_pardiso->nrhs,scale;
581: PetscInt mem = mat_mkl_pardiso->n*mat_mkl_pardiso->nrhs;
583: MatFactorFactorizeSchurComplement(A);
584: /* allocate extra memory if it is needed */
585: scale = 1;
586: if (A->schur_status == MAT_FACTOR_SCHUR_INVERTED) scale = 2;
587: mem *= scale;
588: if (mem > mat_mkl_pardiso->schur_work_size) {
589: o_schur_work = mat_mkl_pardiso->schur_work;
590: PetscMalloc1(mem,&mat_mkl_pardiso->schur_work);
591: }
592: /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
593: if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;
594: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
595: MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
596: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
597: } else { /* if we are solving for the interior problem, any value in barray[schur,n] forward-substituted to xarray[schur,n] will be neglected */
598: PetscInt i,n,m=0;
599: for (n=0;n<mat_mkl_pardiso->nrhs;n++) {
600: for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
601: xarray[mat_mkl_pardiso->schur_idxs[i]+m] = 0.;
602: }
603: m += mat_mkl_pardiso->n;
604: }
605: }
607: /* expansion phase */
608: mat_mkl_pardiso->iparm[6-1] = 1;
609: mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
610: MKL_PARDISO (mat_mkl_pardiso->pt,
611: &mat_mkl_pardiso->maxfct,
612: &mat_mkl_pardiso->mnum,
613: &mat_mkl_pardiso->mtype,
614: &mat_mkl_pardiso->phase,
615: &mat_mkl_pardiso->n,
616: mat_mkl_pardiso->a,
617: mat_mkl_pardiso->ia,
618: mat_mkl_pardiso->ja,
619: mat_mkl_pardiso->perm,
620: &mat_mkl_pardiso->nrhs,
621: mat_mkl_pardiso->iparm,
622: &mat_mkl_pardiso->msglvl,
623: (void*)xarray,
624: (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
625: &mat_mkl_pardiso->err);
626: if (o_schur_work) { /* restore original schur_work (minimal size) */
627: PetscFree(mat_mkl_pardiso->schur_work);
628: mat_mkl_pardiso->schur_work = o_schur_work;
629: }
630: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
631: mat_mkl_pardiso->iparm[6-1] = 0;
632: }
633: MatDenseRestoreArrayWrite(X,&xarray);
634: }
635: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
636: return(0);
637: }
639: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info)
640: {
641: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->data;
642: PetscErrorCode ierr;
645: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
646: (*mat_mkl_pardiso->Convert)(A,mat_mkl_pardiso->needsym,MAT_REUSE_MATRIX,&mat_mkl_pardiso->freeaij,&mat_mkl_pardiso->nz,&mat_mkl_pardiso->ia,&mat_mkl_pardiso->ja,(PetscScalar**)&mat_mkl_pardiso->a);
648: mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
649: MKL_PARDISO (mat_mkl_pardiso->pt,
650: &mat_mkl_pardiso->maxfct,
651: &mat_mkl_pardiso->mnum,
652: &mat_mkl_pardiso->mtype,
653: &mat_mkl_pardiso->phase,
654: &mat_mkl_pardiso->n,
655: mat_mkl_pardiso->a,
656: mat_mkl_pardiso->ia,
657: mat_mkl_pardiso->ja,
658: mat_mkl_pardiso->perm,
659: &mat_mkl_pardiso->nrhs,
660: mat_mkl_pardiso->iparm,
661: &mat_mkl_pardiso->msglvl,
662: NULL,
663: (void*)mat_mkl_pardiso->schur,
664: &mat_mkl_pardiso->err);
665: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
667: /* report flops */
668: if (mat_mkl_pardiso->iparm[18] > 0) {
669: PetscLogFlops(PetscPowRealInt(10.,6)*mat_mkl_pardiso->iparm[18]);
670: }
672: if (F->schur) { /* schur output from pardiso is in row major format */
673: #if defined(PETSC_HAVE_CUDA)
674: F->schur->offloadmask = PETSC_OFFLOAD_CPU;
675: #endif
676: MatFactorRestoreSchurComplement(F,NULL,MAT_FACTOR_SCHUR_UNFACTORED);
677: MatTranspose(F->schur,MAT_INPLACE_MATRIX,&F->schur);
678: }
679: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
680: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
681: return(0);
682: }
684: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A)
685: {
686: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
687: PetscErrorCode ierr;
688: PetscInt icntl,bs,threads=1;
689: PetscBool flg;
692: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");
694: PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use within PARDISO","None",threads,&threads,&flg);
695: if (flg) PetscSetMKL_PARDISOThreads((int)threads);
697: 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);
698: if (flg) mat_mkl_pardiso->maxfct = icntl;
700: PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);
701: if (flg) mat_mkl_pardiso->mnum = icntl;
703: PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);
704: if (flg) mat_mkl_pardiso->msglvl = icntl;
706: PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);
707: if (flg) {
708: void *pt[IPARM_SIZE];
709: mat_mkl_pardiso->mtype = icntl;
710: icntl = mat_mkl_pardiso->iparm[34];
711: bs = mat_mkl_pardiso->iparm[36];
712: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
713: #if defined(PETSC_USE_REAL_SINGLE)
714: mat_mkl_pardiso->iparm[27] = 1;
715: #else
716: mat_mkl_pardiso->iparm[27] = 0;
717: #endif
718: mat_mkl_pardiso->iparm[34] = icntl;
719: mat_mkl_pardiso->iparm[36] = bs;
720: }
722: PetscOptionsInt("-mat_mkl_pardiso_1","Use default values (if 0)","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);
723: if (flg) mat_mkl_pardiso->iparm[0] = icntl;
725: PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);
726: if (flg) mat_mkl_pardiso->iparm[1] = icntl;
728: PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);
729: if (flg) mat_mkl_pardiso->iparm[3] = icntl;
731: PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);
732: if (flg) mat_mkl_pardiso->iparm[4] = icntl;
734: PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);
735: if (flg) mat_mkl_pardiso->iparm[5] = icntl;
737: PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);
738: if (flg) mat_mkl_pardiso->iparm[7] = icntl;
740: PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);
741: if (flg) mat_mkl_pardiso->iparm[9] = icntl;
743: PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);
744: if (flg) mat_mkl_pardiso->iparm[10] = icntl;
746: PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);
747: if (flg) mat_mkl_pardiso->iparm[11] = icntl;
749: PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);
750: if (flg) mat_mkl_pardiso->iparm[12] = icntl;
752: PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);
753: if (flg) mat_mkl_pardiso->iparm[17] = icntl;
755: PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations (0 to disable)","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);
756: if (flg) mat_mkl_pardiso->iparm[18] = icntl;
758: PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);
759: if (flg) mat_mkl_pardiso->iparm[20] = icntl;
761: PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);
762: if (flg) mat_mkl_pardiso->iparm[23] = icntl;
764: PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);
765: if (flg) mat_mkl_pardiso->iparm[24] = icntl;
767: PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);
768: if (flg) mat_mkl_pardiso->iparm[26] = icntl;
770: PetscOptionsInt("-mat_mkl_pardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_pardiso->iparm[30],&icntl,&flg);
771: if (flg) mat_mkl_pardiso->iparm[30] = icntl;
773: PetscOptionsInt("-mat_mkl_pardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_pardiso->iparm[33],&icntl,&flg);
774: if (flg) mat_mkl_pardiso->iparm[33] = icntl;
776: PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
777: if (flg) mat_mkl_pardiso->iparm[59] = icntl;
778: PetscOptionsEnd();
779: return(0);
780: }
782: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
783: {
785: PetscInt i,bs;
786: PetscBool match;
789: for (i=0; i<IPARM_SIZE; i++) mat_mkl_pardiso->iparm[i] = 0;
790: for (i=0; i<IPARM_SIZE; i++) mat_mkl_pardiso->pt[i] = 0;
791: #if defined(PETSC_USE_REAL_SINGLE)
792: mat_mkl_pardiso->iparm[27] = 1;
793: #else
794: mat_mkl_pardiso->iparm[27] = 0;
795: #endif
796: /* Default options for both sym and unsym */
797: mat_mkl_pardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */
798: mat_mkl_pardiso->iparm[ 1] = 2; /* Metis reordering */
799: mat_mkl_pardiso->iparm[ 5] = 0; /* Write solution into x */
800: mat_mkl_pardiso->iparm[ 7] = 0; /* Max number of iterative refinement steps */
801: mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
802: mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
803: #if 0
804: mat_mkl_pardiso->iparm[23] = 1; /* Parallel factorization control*/
805: #endif
806: PetscObjectTypeCompareAny((PetscObject)A,&match,MATSEQBAIJ,MATSEQSBAIJ,"");
807: MatGetBlockSize(A,&bs);
808: if (!match || bs == 1) {
809: mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
810: mat_mkl_pardiso->n = A->rmap->N;
811: } else {
812: mat_mkl_pardiso->iparm[34] = 0; /* Cluster Sparse Solver use Fortran-style indexing for ia and ja arrays */
813: mat_mkl_pardiso->iparm[36] = bs;
814: mat_mkl_pardiso->n = A->rmap->N/bs;
815: }
816: mat_mkl_pardiso->iparm[39] = 0; /* Input: matrix/rhs/solution stored on rank-0 */
818: mat_mkl_pardiso->CleanUp = PETSC_FALSE;
819: mat_mkl_pardiso->maxfct = 1; /* Maximum number of numerical factorizations. */
820: mat_mkl_pardiso->mnum = 1; /* Which factorization to use. */
821: mat_mkl_pardiso->msglvl = 0; /* 0: do not print 1: Print statistical information in file */
822: mat_mkl_pardiso->phase = -1;
823: mat_mkl_pardiso->err = 0;
825: mat_mkl_pardiso->nrhs = 1;
826: mat_mkl_pardiso->err = 0;
827: mat_mkl_pardiso->phase = -1;
829: if (ftype == MAT_FACTOR_LU) {
830: mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
831: mat_mkl_pardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
832: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
833: } else {
834: mat_mkl_pardiso->iparm[ 9] = 8; /* Perturb the pivot elements with 1E-8 */
835: mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
836: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
837: #if defined(PETSC_USE_DEBUG)
838: mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
839: #endif
840: }
841: PetscCalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
842: mat_mkl_pardiso->schur_size = 0;
843: return(0);
844: }
846: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
847: {
848: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
849: PetscErrorCode ierr;
852: mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
853: PetscSetMKL_PARDISOFromOptions(F,A);
854: /* throw away any previously computed structure */
855: if (mat_mkl_pardiso->freeaij) {
856: PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
857: if (mat_mkl_pardiso->iparm[34] == 1) {
858: PetscFree(mat_mkl_pardiso->a);
859: }
860: }
861: (*mat_mkl_pardiso->Convert)(A,mat_mkl_pardiso->needsym,MAT_INITIAL_MATRIX,&mat_mkl_pardiso->freeaij,&mat_mkl_pardiso->nz,&mat_mkl_pardiso->ia,&mat_mkl_pardiso->ja,(PetscScalar**)&mat_mkl_pardiso->a);
862: if (mat_mkl_pardiso->iparm[34] == 1) mat_mkl_pardiso->n = A->rmap->N;
863: else mat_mkl_pardiso->n = A->rmap->N/A->rmap->bs;
865: mat_mkl_pardiso->phase = JOB_ANALYSIS;
867: /* reset flops counting if requested */
868: if (mat_mkl_pardiso->iparm[18]) mat_mkl_pardiso->iparm[18] = -1;
870: MKL_PARDISO (mat_mkl_pardiso->pt,
871: &mat_mkl_pardiso->maxfct,
872: &mat_mkl_pardiso->mnum,
873: &mat_mkl_pardiso->mtype,
874: &mat_mkl_pardiso->phase,
875: &mat_mkl_pardiso->n,
876: mat_mkl_pardiso->a,
877: mat_mkl_pardiso->ia,
878: mat_mkl_pardiso->ja,
879: mat_mkl_pardiso->perm,
880: &mat_mkl_pardiso->nrhs,
881: mat_mkl_pardiso->iparm,
882: &mat_mkl_pardiso->msglvl,
883: NULL,
884: NULL,
885: &mat_mkl_pardiso->err);
886: if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err);
888: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
890: if (F->factortype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
891: else F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO;
893: F->ops->solve = MatSolve_MKL_PARDISO;
894: F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO;
895: F->ops->matsolve = MatMatSolve_MKL_PARDISO;
896: return(0);
897: }
899: PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
900: {
904: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
905: return(0);
906: }
908: #if !defined(PETSC_USE_COMPLEX)
909: PetscErrorCode MatGetInertia_MKL_PARDISO(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
910: {
911: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)F->data;
914: if (nneg) *nneg = mat_mkl_pardiso->iparm[22];
915: if (npos) *npos = mat_mkl_pardiso->iparm[21];
916: if (nzero) *nzero = F->rmap->N - (mat_mkl_pardiso->iparm[22] + mat_mkl_pardiso->iparm[21]);
917: return(0);
918: }
919: #endif
921: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,const MatFactorInfo *info)
922: {
926: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
927: #if defined(PETSC_USE_COMPLEX)
928: F->ops->getinertia = NULL;
929: #else
930: F->ops->getinertia = MatGetInertia_MKL_PARDISO;
931: #endif
932: return(0);
933: }
935: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
936: {
937: PetscErrorCode ierr;
938: PetscBool iascii;
939: PetscViewerFormat format;
940: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
941: PetscInt i;
944: if (A->ops->solve != MatSolve_MKL_PARDISO) return(0);
946: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
947: if (iascii) {
948: PetscViewerGetFormat(viewer,&format);
949: if (format == PETSC_VIEWER_ASCII_INFO) {
950: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
951: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);
952: for (i=1; i<=64; i++) {
953: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
954: }
955: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);
956: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);
957: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);
958: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);
959: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);
960: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);
961: }
962: }
963: return(0);
964: }
966: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
967: {
968: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)A->data;
971: info->block_size = 1.0;
972: info->nz_used = mat_mkl_pardiso->iparm[17];
973: info->nz_allocated = mat_mkl_pardiso->iparm[17];
974: info->nz_unneeded = 0.0;
975: info->assemblies = 0.0;
976: info->mallocs = 0.0;
977: info->memory = 0.0;
978: info->fill_ratio_given = 0;
979: info->fill_ratio_needed = 0;
980: info->factor_mallocs = 0;
981: return(0);
982: }
984: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
985: {
986: PetscInt backup,bs;
987: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
990: if (icntl <= 64) {
991: mat_mkl_pardiso->iparm[icntl - 1] = ival;
992: } else {
993: if (icntl == 65) PetscSetMKL_PARDISOThreads(ival);
994: else if (icntl == 66) mat_mkl_pardiso->maxfct = ival;
995: else if (icntl == 67) mat_mkl_pardiso->mnum = ival;
996: else if (icntl == 68) mat_mkl_pardiso->msglvl = ival;
997: else if (icntl == 69) {
998: void *pt[IPARM_SIZE];
999: backup = mat_mkl_pardiso->iparm[34];
1000: bs = mat_mkl_pardiso->iparm[36];
1001: mat_mkl_pardiso->mtype = ival;
1002: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
1003: #if defined(PETSC_USE_REAL_SINGLE)
1004: mat_mkl_pardiso->iparm[27] = 1;
1005: #else
1006: mat_mkl_pardiso->iparm[27] = 0;
1007: #endif
1008: mat_mkl_pardiso->iparm[34] = backup;
1009: mat_mkl_pardiso->iparm[36] = bs;
1010: } else if (icntl==70) mat_mkl_pardiso->solve_interior = (PetscBool)!!ival;
1011: }
1012: return(0);
1013: }
1015: /*@
1016: MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters
1018: Logically Collective on Mat
1020: Input Parameters:
1021: + F - the factored matrix obtained by calling MatGetFactor()
1022: . icntl - index of Mkl_Pardiso parameter
1023: - ival - value of Mkl_Pardiso parameter
1025: Options Database:
1026: . -mat_mkl_pardiso_<icntl> <ival>
1028: Level: beginner
1030: References:
1031: . Mkl_Pardiso Users' Guide
1033: .seealso: MatGetFactor()
1034: @*/
1035: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
1036: {
1040: PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
1041: return(0);
1042: }
1044: /*MC
1045: MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for
1046: sequential matrices via the external package MKL_PARDISO.
1048: Works with MATSEQAIJ matrices
1050: Use -pc_type lu -pc_factor_mat_solver_type mkl_pardiso to use this direct solver
1052: Options Database Keys:
1053: + -mat_mkl_pardiso_65 - Number of threads to use within MKL_PARDISO
1054: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
1055: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
1056: . -mat_mkl_pardiso_68 - Message level information
1057: . -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
1058: . -mat_mkl_pardiso_1 - Use default values
1059: . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix
1060: . -mat_mkl_pardiso_4 - Preconditioned CGS/CG
1061: . -mat_mkl_pardiso_5 - User permutation
1062: . -mat_mkl_pardiso_6 - Write solution on x
1063: . -mat_mkl_pardiso_8 - Iterative refinement step
1064: . -mat_mkl_pardiso_10 - Pivoting perturbation
1065: . -mat_mkl_pardiso_11 - Scaling vectors
1066: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
1067: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
1068: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
1069: . -mat_mkl_pardiso_19 - Report number of floating point operations
1070: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
1071: . -mat_mkl_pardiso_24 - Parallel factorization control
1072: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
1073: . -mat_mkl_pardiso_27 - Matrix checker
1074: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
1075: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
1076: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode
1078: Level: beginner
1080: For more information please check mkl_pardiso manual
1082: .seealso: PCFactorSetMatSolverType(), MatSolverType
1084: M*/
1085: static PetscErrorCode MatFactorGetSolverType_mkl_pardiso(Mat A, MatSolverType *type)
1086: {
1088: *type = MATSOLVERMKL_PARDISO;
1089: return(0);
1090: }
1092: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
1093: {
1094: Mat B;
1095: PetscErrorCode ierr;
1096: Mat_MKL_PARDISO *mat_mkl_pardiso;
1097: PetscBool isSeqAIJ,isSeqBAIJ,isSeqSBAIJ;
1100: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1101: PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1102: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1103: MatCreate(PetscObjectComm((PetscObject)A),&B);
1104: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1105: PetscStrallocpy("mkl_pardiso",&((PetscObject)B)->type_name);
1106: MatSetUp(B);
1108: PetscNewLog(B,&mat_mkl_pardiso);
1109: B->data = mat_mkl_pardiso;
1111: MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
1112: if (ftype == MAT_FACTOR_LU) {
1113: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
1114: B->factortype = MAT_FACTOR_LU;
1115: mat_mkl_pardiso->needsym = PETSC_FALSE;
1116: if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1117: else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1118: else if (isSeqSBAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU factor with SEQSBAIJ format! Use MAT_FACTOR_CHOLESKY instead");
1119: else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU with %s format",((PetscObject)A)->type_name);
1120: #if defined(PETSC_USE_COMPLEX)
1121: mat_mkl_pardiso->mtype = 13;
1122: #else
1123: mat_mkl_pardiso->mtype = 11;
1124: #endif
1125: } else {
1126: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
1127: B->factortype = MAT_FACTOR_CHOLESKY;
1128: if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1129: else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1130: else if (isSeqSBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqsbaij;
1131: else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with %s format",((PetscObject)A)->type_name);
1133: mat_mkl_pardiso->needsym = PETSC_TRUE;
1134: #if !defined(PETSC_USE_COMPLEX)
1135: if (A->spd_set && A->spd) mat_mkl_pardiso->mtype = 2;
1136: else mat_mkl_pardiso->mtype = -2;
1137: #else
1138: mat_mkl_pardiso->mtype = 6;
1139: if (A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with Hermitian matrices! Use MAT_FACTOR_LU instead");
1140: #endif
1141: }
1142: B->ops->destroy = MatDestroy_MKL_PARDISO;
1143: B->ops->view = MatView_MKL_PARDISO;
1144: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
1145: B->factortype = ftype;
1146: B->assembled = PETSC_TRUE;
1148: PetscFree(B->solvertype);
1149: PetscStrallocpy(MATSOLVERMKL_PARDISO,&B->solvertype);
1151: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_pardiso);
1152: PetscObjectComposeFunction((PetscObject)B,"MatFactorSetSchurIS_C",MatFactorSetSchurIS_MKL_PARDISO);
1153: PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
1155: *F = B;
1156: return(0);
1157: }
1159: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void)
1160: {
1164: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1165: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1166: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQBAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1167: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1168: return(0);
1169: }