Actual source code: mkl_pardiso.c
petsc-3.14.6 2021-03-30
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
57: /*
58: * Internal data structure.
59: * For more information check mkl_pardiso manual.
60: */
61: typedef struct {
63: /* Configuration vector*/
64: INT_TYPE iparm[IPARM_SIZE];
66: /*
67: * Internal mkl_pardiso memory location.
68: * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak.
69: */
70: void *pt[IPARM_SIZE];
72: /* Basic mkl_pardiso info*/
73: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
75: /* Matrix structure*/
76: void *a;
77: INT_TYPE *ia, *ja;
79: /* Number of non-zero elements*/
80: INT_TYPE nz;
82: /* Row permutaton vector*/
83: INT_TYPE *perm;
85: /* Define if matrix preserves sparse structure.*/
86: MatStructure matstruc;
88: PetscBool needsym;
89: PetscBool freeaij;
91: /* Schur complement */
92: PetscScalar *schur;
93: PetscInt schur_size;
94: PetscInt *schur_idxs;
95: PetscScalar *schur_work;
96: PetscBLASInt schur_work_size;
97: PetscBool solve_interior;
99: /* True if mkl_pardiso function have been used.*/
100: PetscBool CleanUp;
102: /* Conversion to a format suitable for MKL */
103: PetscErrorCode (*Convert)(Mat, PetscBool, MatReuse, PetscBool*, INT_TYPE*, INT_TYPE**, INT_TYPE**, PetscScalar**);
104: } Mat_MKL_PARDISO;
106: PetscErrorCode MatMKLPardiso_Convert_seqsbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
107: {
108: Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data;
109: PetscInt bs = A->rmap->bs,i;
113: if (!sym) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
114: *v = aa->a;
115: if (bs == 1) { /* already in the correct format */
116: /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
117: *r = (INT_TYPE*)aa->i;
118: *c = (INT_TYPE*)aa->j;
119: *nnz = (INT_TYPE)aa->nz;
120: *free = PETSC_FALSE;
121: } else if (reuse == MAT_INITIAL_MATRIX) {
122: PetscInt m = A->rmap->n,nz = aa->nz;
123: PetscInt *row,*col;
124: PetscMalloc2(m+1,&row,nz,&col);
125: for (i=0; i<m+1; i++) {
126: row[i] = aa->i[i]+1;
127: }
128: for (i=0; i<nz; i++) {
129: col[i] = aa->j[i]+1;
130: }
131: *r = (INT_TYPE*)row;
132: *c = (INT_TYPE*)col;
133: *nnz = (INT_TYPE)nz;
134: *free = PETSC_TRUE;
135: }
136: return(0);
137: }
139: PetscErrorCode MatMKLPardiso_Convert_seqbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
140: {
141: Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)A->data;
142: PetscInt bs = A->rmap->bs,i;
146: if (!sym) {
147: *v = aa->a;
148: if (bs == 1) { /* already in the correct format */
149: /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
150: *r = (INT_TYPE*)aa->i;
151: *c = (INT_TYPE*)aa->j;
152: *nnz = (INT_TYPE)aa->nz;
153: *free = PETSC_FALSE;
154: return(0);
155: } else if (reuse == MAT_INITIAL_MATRIX) {
156: PetscInt m = A->rmap->n,nz = aa->nz;
157: PetscInt *row,*col;
158: PetscMalloc2(m+1,&row,nz,&col);
159: for (i=0; i<m+1; i++) {
160: row[i] = aa->i[i]+1;
161: }
162: for (i=0; i<nz; i++) {
163: col[i] = aa->j[i]+1;
164: }
165: *r = (INT_TYPE*)row;
166: *c = (INT_TYPE*)col;
167: *nnz = (INT_TYPE)nz;
168: }
169: *free = PETSC_TRUE;
170: } else {
171: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
172: }
173: return(0);
174: }
176: PetscErrorCode MatMKLPardiso_Convert_seqaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v)
177: {
178: Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data;
179: PetscScalar *aav;
183: MatSeqAIJGetArrayRead(A,(const PetscScalar**)&aav);
184: if (!sym) { /* already in the correct format */
185: *v = aav;
186: *r = (INT_TYPE*)aa->i;
187: *c = (INT_TYPE*)aa->j;
188: *nnz = (INT_TYPE)aa->nz;
189: *free = PETSC_FALSE;
190: } else if (reuse == MAT_INITIAL_MATRIX) { /* need to get the triangular part */
191: PetscScalar *vals,*vv;
192: PetscInt *row,*col,*jj;
193: PetscInt m = A->rmap->n,nz,i;
195: nz = 0;
196: for (i=0; i<m; i++) nz += aa->i[i+1] - aa->diag[i];
197: PetscMalloc2(m+1,&row,nz,&col);
198: PetscMalloc1(nz,&vals);
199: jj = col;
200: vv = vals;
202: row[0] = 0;
203: for (i=0; i<m; i++) {
204: PetscInt *aj = aa->j + aa->diag[i];
205: PetscScalar *av = aav + aa->diag[i];
206: PetscInt rl = aa->i[i+1] - aa->diag[i],j;
208: for (j=0; j<rl; j++) {
209: *jj = *aj; jj++; aj++;
210: *vv = *av; vv++; av++;
211: }
212: row[i+1] = row[i] + rl;
213: }
214: *v = vals;
215: *r = (INT_TYPE*)row;
216: *c = (INT_TYPE*)col;
217: *nnz = (INT_TYPE)nz;
218: *free = PETSC_TRUE;
219: } else {
220: PetscScalar *vv;
221: PetscInt m = A->rmap->n,i;
223: vv = *v;
224: for (i=0; i<m; i++) {
225: PetscScalar *av = aav + aa->diag[i];
226: PetscInt rl = aa->i[i+1] - aa->diag[i],j;
227: for (j=0; j<rl; j++) {
228: *vv = *av; vv++; av++;
229: }
230: }
231: *free = PETSC_TRUE;
232: }
233: MatSeqAIJRestoreArrayRead(A,(const PetscScalar**)&aav);
234: return(0);
235: }
238: static PetscErrorCode MatMKLPardisoSolveSchur_Private(Mat F, PetscScalar *B, PetscScalar *X)
239: {
240: Mat_MKL_PARDISO *mpardiso = (Mat_MKL_PARDISO*)F->data;
241: Mat S,Xmat,Bmat;
242: MatFactorSchurStatus schurstatus;
243: PetscErrorCode ierr;
246: MatFactorFactorizeSchurComplement(F);
247: MatFactorGetSchurComplement(F,&S,&schurstatus);
248: if (X == B && schurstatus == MAT_FACTOR_SCHUR_INVERTED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"X and B cannot point to the same address");
249: MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,B,&Bmat);
250: MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,X,&Xmat);
251: MatSetType(Bmat,((PetscObject)S)->type_name);
252: MatSetType(Xmat,((PetscObject)S)->type_name);
253: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
254: MatBindToCPU(Xmat,S->boundtocpu);
255: MatBindToCPU(Bmat,S->boundtocpu);
256: #endif
258: #if defined(PETSC_USE_COMPLEX)
259: if (mpardiso->iparm[12-1] == 1) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Hermitian solve not implemented yet");
260: #endif
262: switch (schurstatus) {
263: case MAT_FACTOR_SCHUR_FACTORED:
264: if (!mpardiso->iparm[12-1]) {
265: MatMatSolve(S,Bmat,Xmat);
266: } else { /* transpose solve */
267: MatMatSolveTranspose(S,Bmat,Xmat);
268: }
269: break;
270: case MAT_FACTOR_SCHUR_INVERTED:
271: MatProductCreateWithMat(S,Bmat,NULL,Xmat);
272: if (!mpardiso->iparm[12-1]) {
273: MatProductSetType(Xmat,MATPRODUCT_AB);
274: } else { /* transpose solve */
275: MatProductSetType(Xmat,MATPRODUCT_AtB);
276: }
277: MatProductSetFromOptions(Xmat);
278: MatProductSymbolic(Xmat);
279: MatProductNumeric(Xmat);
280: MatProductClear(Xmat);
281: break;
282: default:
283: SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status);
284: break;
285: }
286: MatFactorRestoreSchurComplement(F,&S,schurstatus);
287: MatDestroy(&Bmat);
288: MatDestroy(&Xmat);
289: return(0);
290: }
292: PetscErrorCode MatFactorSetSchurIS_MKL_PARDISO(Mat F, IS is)
293: {
294: Mat_MKL_PARDISO *mpardiso = (Mat_MKL_PARDISO*)F->data;
295: const PetscScalar *arr;
296: const PetscInt *idxs;
297: PetscInt size,i;
298: PetscMPIInt csize;
299: PetscBool sorted;
300: PetscErrorCode ierr;
303: MPI_Comm_size(PetscObjectComm((PetscObject)F),&csize);
304: if (csize > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MKL_PARDISO parallel Schur complements not yet supported from PETSc");
305: ISSorted(is,&sorted);
306: if (!sorted) {
307: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS for MKL_PARDISO Schur complements needs to be sorted");
308: }
309: ISGetLocalSize(is,&size);
310: PetscFree(mpardiso->schur_work);
311: PetscBLASIntCast(PetscMax(mpardiso->n,2*size),&mpardiso->schur_work_size);
312: PetscMalloc1(mpardiso->schur_work_size,&mpardiso->schur_work);
313: MatDestroy(&F->schur);
314: MatCreateSeqDense(PETSC_COMM_SELF,size,size,NULL,&F->schur);
315: MatDenseGetArrayRead(F->schur,&arr);
316: mpardiso->schur = (PetscScalar*)arr;
317: mpardiso->schur_size = size;
318: MatDenseRestoreArrayRead(F->schur,&arr);
319: if (mpardiso->mtype == 2) {
320: MatSetOption(F->schur,MAT_SPD,PETSC_TRUE);
321: }
323: PetscFree(mpardiso->schur_idxs);
324: PetscMalloc1(size,&mpardiso->schur_idxs);
325: PetscArrayzero(mpardiso->perm,mpardiso->n);
326: ISGetIndices(is,&idxs);
327: PetscArraycpy(mpardiso->schur_idxs,idxs,size);
328: for (i=0;i<size;i++) mpardiso->perm[idxs[i]] = 1;
329: ISRestoreIndices(is,&idxs);
330: if (size) { /* turn on Schur switch if the set of indices is not empty */
331: mpardiso->iparm[36-1] = 2;
332: }
333: return(0);
334: }
336: PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
337: {
338: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
339: PetscErrorCode ierr;
342: if (mat_mkl_pardiso->CleanUp) {
343: mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
345: MKL_PARDISO (mat_mkl_pardiso->pt,
346: &mat_mkl_pardiso->maxfct,
347: &mat_mkl_pardiso->mnum,
348: &mat_mkl_pardiso->mtype,
349: &mat_mkl_pardiso->phase,
350: &mat_mkl_pardiso->n,
351: NULL,
352: NULL,
353: NULL,
354: NULL,
355: &mat_mkl_pardiso->nrhs,
356: mat_mkl_pardiso->iparm,
357: &mat_mkl_pardiso->msglvl,
358: NULL,
359: NULL,
360: &mat_mkl_pardiso->err);
361: }
362: PetscFree(mat_mkl_pardiso->perm);
363: PetscFree(mat_mkl_pardiso->schur_work);
364: PetscFree(mat_mkl_pardiso->schur_idxs);
365: if (mat_mkl_pardiso->freeaij) {
366: PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
367: if (mat_mkl_pardiso->iparm[34] == 1) {
368: PetscFree(mat_mkl_pardiso->a);
369: }
370: }
371: PetscFree(A->data);
373: /* clear composed functions */
374: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
375: PetscObjectComposeFunction((PetscObject)A,"MatFactorSetSchurIS_C",NULL);
376: PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);
377: return(0);
378: }
380: static PetscErrorCode MatMKLPardisoScatterSchur_Private(Mat_MKL_PARDISO *mpardiso, PetscScalar *whole, PetscScalar *schur, PetscBool reduce)
381: {
383: if (reduce) { /* data given for the whole matrix */
384: PetscInt i,m=0,p=0;
385: for (i=0;i<mpardiso->nrhs;i++) {
386: PetscInt j;
387: for (j=0;j<mpardiso->schur_size;j++) {
388: schur[p+j] = whole[m+mpardiso->schur_idxs[j]];
389: }
390: m += mpardiso->n;
391: p += mpardiso->schur_size;
392: }
393: } else { /* from Schur to whole */
394: PetscInt i,m=0,p=0;
395: for (i=0;i<mpardiso->nrhs;i++) {
396: PetscInt j;
397: for (j=0;j<mpardiso->schur_size;j++) {
398: whole[m+mpardiso->schur_idxs[j]] = schur[p+j];
399: }
400: m += mpardiso->n;
401: p += mpardiso->schur_size;
402: }
403: }
404: return(0);
405: }
407: PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x)
408: {
409: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
410: PetscErrorCode ierr;
411: PetscScalar *xarray;
412: const PetscScalar *barray;
415: mat_mkl_pardiso->nrhs = 1;
416: VecGetArray(x,&xarray);
417: VecGetArrayRead(b,&barray);
419: if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
420: else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
422: if (barray == xarray) { /* if the two vectors share the same memory */
423: PetscScalar *work;
424: if (!mat_mkl_pardiso->schur_work) {
425: PetscMalloc1(mat_mkl_pardiso->n,&work);
426: } else {
427: work = mat_mkl_pardiso->schur_work;
428: }
429: mat_mkl_pardiso->iparm[6-1] = 1;
430: MKL_PARDISO (mat_mkl_pardiso->pt,
431: &mat_mkl_pardiso->maxfct,
432: &mat_mkl_pardiso->mnum,
433: &mat_mkl_pardiso->mtype,
434: &mat_mkl_pardiso->phase,
435: &mat_mkl_pardiso->n,
436: mat_mkl_pardiso->a,
437: mat_mkl_pardiso->ia,
438: mat_mkl_pardiso->ja,
439: NULL,
440: &mat_mkl_pardiso->nrhs,
441: mat_mkl_pardiso->iparm,
442: &mat_mkl_pardiso->msglvl,
443: (void*)xarray,
444: (void*)work,
445: &mat_mkl_pardiso->err);
446: if (!mat_mkl_pardiso->schur_work) {
447: PetscFree(work);
448: }
449: } else {
450: mat_mkl_pardiso->iparm[6-1] = 0;
451: MKL_PARDISO (mat_mkl_pardiso->pt,
452: &mat_mkl_pardiso->maxfct,
453: &mat_mkl_pardiso->mnum,
454: &mat_mkl_pardiso->mtype,
455: &mat_mkl_pardiso->phase,
456: &mat_mkl_pardiso->n,
457: mat_mkl_pardiso->a,
458: mat_mkl_pardiso->ia,
459: mat_mkl_pardiso->ja,
460: mat_mkl_pardiso->perm,
461: &mat_mkl_pardiso->nrhs,
462: mat_mkl_pardiso->iparm,
463: &mat_mkl_pardiso->msglvl,
464: (void*)barray,
465: (void*)xarray,
466: &mat_mkl_pardiso->err);
467: }
468: VecRestoreArrayRead(b,&barray);
470: 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);
472: if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
473: PetscInt shift = mat_mkl_pardiso->schur_size;
475: /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
476: if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;
478: if (!mat_mkl_pardiso->solve_interior) {
479: /* solve Schur complement */
480: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
481: MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
482: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
483: } else { /* if we are solving for the interior problem, any value in barray[schur] forward-substituted to xarray[schur] will be neglected */
484: PetscInt i;
485: for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
486: xarray[mat_mkl_pardiso->schur_idxs[i]] = 0.;
487: }
488: }
490: /* expansion phase */
491: mat_mkl_pardiso->iparm[6-1] = 1;
492: mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
493: MKL_PARDISO (mat_mkl_pardiso->pt,
494: &mat_mkl_pardiso->maxfct,
495: &mat_mkl_pardiso->mnum,
496: &mat_mkl_pardiso->mtype,
497: &mat_mkl_pardiso->phase,
498: &mat_mkl_pardiso->n,
499: mat_mkl_pardiso->a,
500: mat_mkl_pardiso->ia,
501: mat_mkl_pardiso->ja,
502: mat_mkl_pardiso->perm,
503: &mat_mkl_pardiso->nrhs,
504: mat_mkl_pardiso->iparm,
505: &mat_mkl_pardiso->msglvl,
506: (void*)xarray,
507: (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
508: &mat_mkl_pardiso->err);
510: 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);
511: mat_mkl_pardiso->iparm[6-1] = 0;
512: }
513: VecRestoreArray(x,&xarray);
514: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
515: return(0);
516: }
518: PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x)
519: {
520: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
521: PetscInt oiparm12;
522: PetscErrorCode ierr;
525: oiparm12 = mat_mkl_pardiso->iparm[12 - 1];
526: mat_mkl_pardiso->iparm[12 - 1] = 2;
527: MatSolve_MKL_PARDISO(A,b,x);
528: mat_mkl_pardiso->iparm[12 - 1] = oiparm12;
529: return(0);
530: }
532: PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X)
533: {
534: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->data;
535: PetscErrorCode ierr;
536: const PetscScalar *barray;
537: PetscScalar *xarray;
538: PetscBool flg;
541: PetscObjectBaseTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
542: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix");
543: if (X != B) {
544: PetscObjectBaseTypeCompare((PetscObject)X,MATSEQDENSE,&flg);
545: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix");
546: }
548: MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);
550: if (mat_mkl_pardiso->nrhs > 0) {
551: MatDenseGetArrayRead(B,&barray);
552: MatDenseGetArray(X,&xarray);
554: if (barray == xarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"B and X cannot share the same memory location");
555: if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
556: else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
558: MKL_PARDISO (mat_mkl_pardiso->pt,
559: &mat_mkl_pardiso->maxfct,
560: &mat_mkl_pardiso->mnum,
561: &mat_mkl_pardiso->mtype,
562: &mat_mkl_pardiso->phase,
563: &mat_mkl_pardiso->n,
564: mat_mkl_pardiso->a,
565: mat_mkl_pardiso->ia,
566: mat_mkl_pardiso->ja,
567: mat_mkl_pardiso->perm,
568: &mat_mkl_pardiso->nrhs,
569: mat_mkl_pardiso->iparm,
570: &mat_mkl_pardiso->msglvl,
571: (void*)barray,
572: (void*)xarray,
573: &mat_mkl_pardiso->err);
574: 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);
576: MatDenseRestoreArrayRead(B,&barray);
577: if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
578: PetscScalar *o_schur_work = NULL;
579: PetscInt shift = mat_mkl_pardiso->schur_size*mat_mkl_pardiso->nrhs,scale;
580: PetscInt mem = mat_mkl_pardiso->n*mat_mkl_pardiso->nrhs;
582: /* allocate extra memory if it is needed */
583: scale = 1;
584: if (A->schur_status == MAT_FACTOR_SCHUR_INVERTED) scale = 2;
586: mem *= scale;
587: if (mem > mat_mkl_pardiso->schur_work_size) {
588: o_schur_work = mat_mkl_pardiso->schur_work;
589: PetscMalloc1(mem,&mat_mkl_pardiso->schur_work);
590: }
592: /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
593: if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;
595: /* solve Schur complement */
596: if (!mat_mkl_pardiso->solve_interior) {
597: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);
598: MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);
599: MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);
600: } else { /* if we are solving for the interior problem, any value in barray[schur,n] forward-substituted to xarray[schur,n] will be neglected */
601: PetscInt i,n,m=0;
602: for (n=0;n<mat_mkl_pardiso->nrhs;n++) {
603: for (i=0;i<mat_mkl_pardiso->schur_size;i++) {
604: xarray[mat_mkl_pardiso->schur_idxs[i]+m] = 0.;
605: }
606: m += mat_mkl_pardiso->n;
607: }
608: }
610: /* expansion phase */
611: mat_mkl_pardiso->iparm[6-1] = 1;
612: mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
613: MKL_PARDISO (mat_mkl_pardiso->pt,
614: &mat_mkl_pardiso->maxfct,
615: &mat_mkl_pardiso->mnum,
616: &mat_mkl_pardiso->mtype,
617: &mat_mkl_pardiso->phase,
618: &mat_mkl_pardiso->n,
619: mat_mkl_pardiso->a,
620: mat_mkl_pardiso->ia,
621: mat_mkl_pardiso->ja,
622: mat_mkl_pardiso->perm,
623: &mat_mkl_pardiso->nrhs,
624: mat_mkl_pardiso->iparm,
625: &mat_mkl_pardiso->msglvl,
626: (void*)xarray,
627: (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
628: &mat_mkl_pardiso->err);
629: if (o_schur_work) { /* restore original schur_work (minimal size) */
630: PetscFree(mat_mkl_pardiso->schur_work);
631: mat_mkl_pardiso->schur_work = o_schur_work;
632: }
633: 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);
634: mat_mkl_pardiso->iparm[6-1] = 0;
635: }
636: }
637: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
638: return(0);
639: }
641: PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info)
642: {
643: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->data;
644: PetscErrorCode ierr;
647: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
648: (*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);
650: mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
651: MKL_PARDISO (mat_mkl_pardiso->pt,
652: &mat_mkl_pardiso->maxfct,
653: &mat_mkl_pardiso->mnum,
654: &mat_mkl_pardiso->mtype,
655: &mat_mkl_pardiso->phase,
656: &mat_mkl_pardiso->n,
657: mat_mkl_pardiso->a,
658: mat_mkl_pardiso->ia,
659: mat_mkl_pardiso->ja,
660: mat_mkl_pardiso->perm,
661: &mat_mkl_pardiso->nrhs,
662: mat_mkl_pardiso->iparm,
663: &mat_mkl_pardiso->msglvl,
664: NULL,
665: (void*)mat_mkl_pardiso->schur,
666: &mat_mkl_pardiso->err);
667: 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);
669: /* report flops */
670: if (mat_mkl_pardiso->iparm[18] > 0) {
671: PetscLogFlops(PetscPowRealInt(10.,6)*mat_mkl_pardiso->iparm[18]);
672: }
674: if (F->schur) { /* schur output from pardiso is in row major format */
675: #if defined(PETSC_HAVE_CUDA)
676: F->schur->offloadmask = PETSC_OFFLOAD_CPU;
677: #endif
678: MatFactorRestoreSchurComplement(F,NULL,MAT_FACTOR_SCHUR_UNFACTORED);
679: MatTranspose(F->schur,MAT_INPLACE_MATRIX,&F->schur);
680: }
681: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
682: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
683: return(0);
684: }
686: PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A)
687: {
688: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
689: PetscErrorCode ierr;
690: PetscInt icntl,bs,threads=1;
691: PetscBool flg;
694: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");
696: PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use within PARDISO","None",threads,&threads,&flg);
697: if (flg) PetscSetMKL_PARDISOThreads((int)threads);
699: 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);
700: if (flg) mat_mkl_pardiso->maxfct = icntl;
702: PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);
703: if (flg) mat_mkl_pardiso->mnum = icntl;
705: PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);
706: if (flg) mat_mkl_pardiso->msglvl = icntl;
708: PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);
709: if (flg) {
710: void *pt[IPARM_SIZE];
711: mat_mkl_pardiso->mtype = icntl;
712: icntl = mat_mkl_pardiso->iparm[34];
713: bs = mat_mkl_pardiso->iparm[36];
714: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
715: #if defined(PETSC_USE_REAL_SINGLE)
716: mat_mkl_pardiso->iparm[27] = 1;
717: #else
718: mat_mkl_pardiso->iparm[27] = 0;
719: #endif
720: mat_mkl_pardiso->iparm[34] = icntl;
721: mat_mkl_pardiso->iparm[36] = bs;
722: }
723: PetscOptionsInt("-mat_mkl_pardiso_1","Use default values (if 0)","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);
725: if (flg && icntl != 0) {
726: PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);
727: if (flg) mat_mkl_pardiso->iparm[1] = icntl;
729: PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);
730: if (flg) mat_mkl_pardiso->iparm[3] = icntl;
732: PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);
733: if (flg) mat_mkl_pardiso->iparm[4] = icntl;
735: PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);
736: if (flg) mat_mkl_pardiso->iparm[5] = icntl;
738: PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);
739: if (flg) mat_mkl_pardiso->iparm[7] = icntl;
741: PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);
742: if (flg) mat_mkl_pardiso->iparm[9] = icntl;
744: PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);
745: if (flg) mat_mkl_pardiso->iparm[10] = icntl;
747: PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);
748: if (flg) mat_mkl_pardiso->iparm[11] = icntl;
750: PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);
751: if (flg) mat_mkl_pardiso->iparm[12] = icntl;
753: PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);
754: if (flg) mat_mkl_pardiso->iparm[17] = icntl;
756: PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations (0 to disable)","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);
757: if (flg) mat_mkl_pardiso->iparm[18] = icntl;
759: PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);
760: if (flg) mat_mkl_pardiso->iparm[20] = icntl;
762: PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);
763: if (flg) mat_mkl_pardiso->iparm[23] = icntl;
765: PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);
766: if (flg) mat_mkl_pardiso->iparm[24] = icntl;
768: PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);
769: if (flg) mat_mkl_pardiso->iparm[26] = icntl;
771: PetscOptionsInt("-mat_mkl_pardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_pardiso->iparm[30],&icntl,&flg);
772: if (flg) mat_mkl_pardiso->iparm[30] = icntl;
774: PetscOptionsInt("-mat_mkl_pardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_pardiso->iparm[33],&icntl,&flg);
775: if (flg) mat_mkl_pardiso->iparm[33] = icntl;
777: PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);
778: if (flg) mat_mkl_pardiso->iparm[59] = icntl;
779: }
780: PetscOptionsEnd();
781: return(0);
782: }
784: PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
785: {
787: PetscInt i,bs;
788: PetscBool match;
791: for (i=0; i<IPARM_SIZE; i++) mat_mkl_pardiso->iparm[i] = 0;
792: for (i=0; i<IPARM_SIZE; i++) mat_mkl_pardiso->pt[i] = 0;
793: #if defined(PETSC_USE_REAL_SINGLE)
794: mat_mkl_pardiso->iparm[27] = 1;
795: #else
796: mat_mkl_pardiso->iparm[27] = 0;
797: #endif
798: /* Default options for both sym and unsym */
799: mat_mkl_pardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */
800: mat_mkl_pardiso->iparm[ 1] = 2; /* Metis reordering */
801: mat_mkl_pardiso->iparm[ 5] = 0; /* Write solution into x */
802: mat_mkl_pardiso->iparm[ 7] = 0; /* Max number of iterative refinement steps */
803: mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
804: mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
805: #if 0
806: mat_mkl_pardiso->iparm[23] = 1; /* Parallel factorization control*/
807: #endif
808: PetscObjectTypeCompareAny((PetscObject)A,&match,MATSEQBAIJ,MATSEQSBAIJ,"");
809: MatGetBlockSize(A,&bs);
810: if (!match || bs == 1) {
811: mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
812: mat_mkl_pardiso->n = A->rmap->N;
813: } else {
814: mat_mkl_pardiso->iparm[34] = 0; /* Cluster Sparse Solver use Fortran-style indexing for ia and ja arrays */
815: mat_mkl_pardiso->iparm[36] = bs;
816: mat_mkl_pardiso->n = A->rmap->N/bs;
817: }
818: mat_mkl_pardiso->iparm[39] = 0; /* Input: matrix/rhs/solution stored on master */
820: mat_mkl_pardiso->CleanUp = PETSC_FALSE;
821: mat_mkl_pardiso->maxfct = 1; /* Maximum number of numerical factorizations. */
822: mat_mkl_pardiso->mnum = 1; /* Which factorization to use. */
823: mat_mkl_pardiso->msglvl = 0; /* 0: do not print 1: Print statistical information in file */
824: mat_mkl_pardiso->phase = -1;
825: mat_mkl_pardiso->err = 0;
827: mat_mkl_pardiso->nrhs = 1;
828: mat_mkl_pardiso->err = 0;
829: mat_mkl_pardiso->phase = -1;
831: if (ftype == MAT_FACTOR_LU) {
832: mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
833: mat_mkl_pardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
834: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
835: } else {
836: mat_mkl_pardiso->iparm[ 9] = 8; /* Perturb the pivot elements with 1E-8 */
837: mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
838: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
839: #if defined(PETSC_USE_DEBUG)
840: mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
841: #endif
842: }
843: PetscMalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);
844: for (i=0; i<A->rmap->N; i++) {
845: mat_mkl_pardiso->perm[i] = 0;
846: }
847: mat_mkl_pardiso->schur_size = 0;
848: return(0);
849: }
851: PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info)
852: {
853: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
854: PetscErrorCode ierr;
857: mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
858: PetscSetMKL_PARDISOFromOptions(F,A);
860: /* throw away any previously computed structure */
861: if (mat_mkl_pardiso->freeaij) {
862: PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);
863: if (mat_mkl_pardiso->iparm[34] == 1) {
864: PetscFree(mat_mkl_pardiso->a);
865: }
866: }
867: (*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);
868: if (mat_mkl_pardiso->iparm[34] == 1) mat_mkl_pardiso->n = A->rmap->N;
869: else mat_mkl_pardiso->n = A->rmap->N/A->rmap->bs;
871: mat_mkl_pardiso->phase = JOB_ANALYSIS;
873: /* reset flops counting if requested */
874: if (mat_mkl_pardiso->iparm[18]) mat_mkl_pardiso->iparm[18] = -1;
876: MKL_PARDISO (mat_mkl_pardiso->pt,
877: &mat_mkl_pardiso->maxfct,
878: &mat_mkl_pardiso->mnum,
879: &mat_mkl_pardiso->mtype,
880: &mat_mkl_pardiso->phase,
881: &mat_mkl_pardiso->n,
882: mat_mkl_pardiso->a,
883: mat_mkl_pardiso->ia,
884: mat_mkl_pardiso->ja,
885: mat_mkl_pardiso->perm,
886: &mat_mkl_pardiso->nrhs,
887: mat_mkl_pardiso->iparm,
888: &mat_mkl_pardiso->msglvl,
889: NULL,
890: NULL,
891: &mat_mkl_pardiso->err);
892: 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);
894: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
896: if (F->factortype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
897: else F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO;
899: F->ops->solve = MatSolve_MKL_PARDISO;
900: F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO;
901: F->ops->matsolve = MatMatSolve_MKL_PARDISO;
902: return(0);
903: }
905: PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
906: {
910: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
911: return(0);
912: }
914: #if !defined(PETSC_USE_COMPLEX)
915: PetscErrorCode MatGetInertia_MKL_PARDISO(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
916: {
917: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)F->data;
920: if (nneg) *nneg = mat_mkl_pardiso->iparm[22];
921: if (npos) *npos = mat_mkl_pardiso->iparm[21];
922: if (nzero) *nzero = F->rmap->N - (mat_mkl_pardiso->iparm[22] + mat_mkl_pardiso->iparm[21]);
923: return(0);
924: }
925: #endif
927: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,const MatFactorInfo *info)
928: {
932: MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);
933: #if defined(PETSC_USE_COMPLEX)
934: F->ops->getinertia = NULL;
935: #else
936: F->ops->getinertia = MatGetInertia_MKL_PARDISO;
937: #endif
938: return(0);
939: }
941: PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
942: {
943: PetscErrorCode ierr;
944: PetscBool iascii;
945: PetscViewerFormat format;
946: Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data;
947: PetscInt i;
950: if (A->ops->solve != MatSolve_MKL_PARDISO) return(0);
952: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
953: if (iascii) {
954: PetscViewerGetFormat(viewer,&format);
955: if (format == PETSC_VIEWER_ASCII_INFO) {
956: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");
957: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);
958: for (i=1; i<=64; i++) {
959: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);
960: }
961: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);
962: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);
963: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);
964: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);
965: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);
966: PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);
967: }
968: }
969: return(0);
970: }
973: PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
974: {
975: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)A->data;
978: info->block_size = 1.0;
979: info->nz_used = mat_mkl_pardiso->iparm[17];
980: info->nz_allocated = mat_mkl_pardiso->iparm[17];
981: info->nz_unneeded = 0.0;
982: info->assemblies = 0.0;
983: info->mallocs = 0.0;
984: info->memory = 0.0;
985: info->fill_ratio_given = 0;
986: info->fill_ratio_needed = 0;
987: info->factor_mallocs = 0;
988: return(0);
989: }
991: PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival)
992: {
993: PetscInt backup,bs;
994: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data;
997: if (icntl <= 64) {
998: mat_mkl_pardiso->iparm[icntl - 1] = ival;
999: } else {
1000: if (icntl == 65) PetscSetMKL_PARDISOThreads(ival);
1001: else if (icntl == 66) mat_mkl_pardiso->maxfct = ival;
1002: else if (icntl == 67) mat_mkl_pardiso->mnum = ival;
1003: else if (icntl == 68) mat_mkl_pardiso->msglvl = ival;
1004: else if (icntl == 69) {
1005: void *pt[IPARM_SIZE];
1006: backup = mat_mkl_pardiso->iparm[34];
1007: bs = mat_mkl_pardiso->iparm[36];
1008: mat_mkl_pardiso->mtype = ival;
1009: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
1010: #if defined(PETSC_USE_REAL_SINGLE)
1011: mat_mkl_pardiso->iparm[27] = 1;
1012: #else
1013: mat_mkl_pardiso->iparm[27] = 0;
1014: #endif
1015: mat_mkl_pardiso->iparm[34] = backup;
1016: mat_mkl_pardiso->iparm[36] = bs;
1017: } else if (icntl==70) mat_mkl_pardiso->solve_interior = (PetscBool)!!ival;
1018: }
1019: return(0);
1020: }
1022: /*@
1023: MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters
1025: Logically Collective on Mat
1027: Input Parameters:
1028: + F - the factored matrix obtained by calling MatGetFactor()
1029: . icntl - index of Mkl_Pardiso parameter
1030: - ival - value of Mkl_Pardiso parameter
1032: Options Database:
1033: . -mat_mkl_pardiso_<icntl> <ival>
1035: Level: beginner
1037: References:
1038: . Mkl_Pardiso Users' Guide
1040: .seealso: MatGetFactor()
1041: @*/
1042: PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
1043: {
1047: PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
1048: return(0);
1049: }
1051: /*MC
1052: MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for
1053: sequential matrices via the external package MKL_PARDISO.
1055: Works with MATSEQAIJ matrices
1057: Use -pc_type lu -pc_factor_mat_solver_type mkl_pardiso to use this direct solver
1059: Options Database Keys:
1060: + -mat_mkl_pardiso_65 - Number of threads to use within MKL_PARDISO
1061: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
1062: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
1063: . -mat_mkl_pardiso_68 - Message level information
1064: . -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
1065: . -mat_mkl_pardiso_1 - Use default values
1066: . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix
1067: . -mat_mkl_pardiso_4 - Preconditioned CGS/CG
1068: . -mat_mkl_pardiso_5 - User permutation
1069: . -mat_mkl_pardiso_6 - Write solution on x
1070: . -mat_mkl_pardiso_8 - Iterative refinement step
1071: . -mat_mkl_pardiso_10 - Pivoting perturbation
1072: . -mat_mkl_pardiso_11 - Scaling vectors
1073: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
1074: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
1075: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
1076: . -mat_mkl_pardiso_19 - Report number of floating point operations
1077: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
1078: . -mat_mkl_pardiso_24 - Parallel factorization control
1079: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
1080: . -mat_mkl_pardiso_27 - Matrix checker
1081: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
1082: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
1083: - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode
1085: Level: beginner
1087: For more information please check mkl_pardiso manual
1089: .seealso: PCFactorSetMatSolverType(), MatSolverType
1091: M*/
1092: static PetscErrorCode MatFactorGetSolverType_mkl_pardiso(Mat A, MatSolverType *type)
1093: {
1095: *type = MATSOLVERMKL_PARDISO;
1096: return(0);
1097: }
1099: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F)
1100: {
1101: Mat B;
1102: PetscErrorCode ierr;
1103: Mat_MKL_PARDISO *mat_mkl_pardiso;
1104: PetscBool isSeqAIJ,isSeqBAIJ,isSeqSBAIJ;
1107: PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
1108: PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);
1109: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);
1110: MatCreate(PetscObjectComm((PetscObject)A),&B);
1111: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
1112: PetscStrallocpy("mkl_pardiso",&((PetscObject)B)->type_name);
1113: MatSetUp(B);
1115: PetscNewLog(B,&mat_mkl_pardiso);
1116: B->data = mat_mkl_pardiso;
1118: MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);
1119: if (ftype == MAT_FACTOR_LU) {
1120: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
1121: B->factortype = MAT_FACTOR_LU;
1122: mat_mkl_pardiso->needsym = PETSC_FALSE;
1123: if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1124: else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1125: else if (isSeqSBAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU factor with SEQSBAIJ format! Use MAT_FACTOR_CHOLESKY instead");
1126: else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU with %s format",((PetscObject)A)->type_name);
1127: #if defined(PETSC_USE_COMPLEX)
1128: mat_mkl_pardiso->mtype = 13;
1129: #else
1130: mat_mkl_pardiso->mtype = 11;
1131: #endif
1132: } else {
1133: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
1134: B->factortype = MAT_FACTOR_CHOLESKY;
1135: if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
1136: else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
1137: else if (isSeqSBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqsbaij;
1138: else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with %s format",((PetscObject)A)->type_name);
1140: mat_mkl_pardiso->needsym = PETSC_TRUE;
1141: #if !defined(PETSC_USE_COMPLEX)
1142: if (A->spd_set && A->spd) mat_mkl_pardiso->mtype = 2;
1143: else mat_mkl_pardiso->mtype = -2;
1144: #else
1145: mat_mkl_pardiso->mtype = 6;
1146: if (A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with Hermitian matrices! Use MAT_FACTOR_LU instead");
1147: #endif
1148: }
1149: B->ops->destroy = MatDestroy_MKL_PARDISO;
1150: B->ops->view = MatView_MKL_PARDISO;
1151: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
1152: B->factortype = ftype;
1153: B->assembled = PETSC_TRUE;
1155: PetscFree(B->solvertype);
1156: PetscStrallocpy(MATSOLVERMKL_PARDISO,&B->solvertype);
1158: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_pardiso);
1159: PetscObjectComposeFunction((PetscObject)B,"MatFactorSetSchurIS_C",MatFactorSetSchurIS_MKL_PARDISO);
1160: PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);
1162: *F = B;
1163: return(0);
1164: }
1166: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void)
1167: {
1171: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1172: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1173: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQBAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);
1174: MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);
1175: return(0);
1176: }