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: */
59: typedef struct {
60: /* Configuration vector*/
61: INT_TYPE iparm[IPARM_SIZE];
63: /*
64: Internal MKL PARDISO memory location.
65: After the first call to MKL PARDISO do not modify pt, as that could cause a serious memory leak.
66: */
67: void *pt[IPARM_SIZE];
69: /* Basic MKL PARDISO info */
70: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
72: /* Matrix structure*/
73: void *a;
74: INT_TYPE *ia, *ja;
76: /* Number of non-zero elements*/
77: INT_TYPE nz;
79: /* Row permutaton vector*/
80: INT_TYPE *perm;
82: /* Define if matrix preserves sparse structure.*/
83: MatStructure matstruc;
85: PetscBool needsym;
86: PetscBool freeaij;
88: /* Schur complement */
89: PetscScalar *schur;
90: PetscInt schur_size;
91: PetscInt *schur_idxs;
92: PetscScalar *schur_work;
93: PetscBLASInt schur_work_size;
94: PetscBool solve_interior;
96: /* True if MKL PARDISO function have been used. */
97: PetscBool CleanUp;
99: /* Conversion to a format suitable for MKL */
100: PetscErrorCode (*Convert)(Mat, PetscBool, MatReuse, PetscBool *, INT_TYPE *, INT_TYPE **, INT_TYPE **, PetscScalar **);
101: } Mat_MKL_PARDISO;
103: static PetscErrorCode MatMKLPardiso_Convert_seqsbaij(Mat A, PetscBool sym, MatReuse reuse, PetscBool *free, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, PetscScalar **v)
104: {
105: Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ *)A->data;
106: PetscInt bs = A->rmap->bs, i;
108: PetscFunctionBegin;
109: PetscCheck(sym, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen");
110: *v = aa->a;
111: if (bs == 1) { /* already in the correct format */
112: /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
113: *r = (INT_TYPE *)aa->i;
114: *c = (INT_TYPE *)aa->j;
115: *nnz = (INT_TYPE)aa->nz;
116: *free = PETSC_FALSE;
117: } else if (reuse == MAT_INITIAL_MATRIX) {
118: PetscInt m = A->rmap->n, nz = aa->nz;
119: PetscInt *row, *col;
120: PetscCall(PetscMalloc2(m + 1, &row, nz, &col));
121: for (i = 0; i < m + 1; i++) row[i] = aa->i[i] + 1;
122: for (i = 0; i < nz; i++) col[i] = aa->j[i] + 1;
123: *r = (INT_TYPE *)row;
124: *c = (INT_TYPE *)col;
125: *nnz = (INT_TYPE)nz;
126: *free = PETSC_TRUE;
127: }
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode MatMKLPardiso_Convert_seqbaij(Mat A, PetscBool sym, MatReuse reuse, PetscBool *free, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, PetscScalar **v)
132: {
133: Mat_SeqBAIJ *aa = (Mat_SeqBAIJ *)A->data;
134: PetscInt bs = A->rmap->bs, i;
136: PetscFunctionBegin;
137: if (!sym) {
138: *v = aa->a;
139: if (bs == 1) { /* already in the correct format */
140: /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */
141: *r = (INT_TYPE *)aa->i;
142: *c = (INT_TYPE *)aa->j;
143: *nnz = (INT_TYPE)aa->nz;
144: *free = PETSC_FALSE;
145: PetscFunctionReturn(PETSC_SUCCESS);
146: } else if (reuse == MAT_INITIAL_MATRIX) {
147: PetscInt m = A->rmap->n, nz = aa->nz;
148: PetscInt *row, *col;
149: PetscCall(PetscMalloc2(m + 1, &row, nz, &col));
150: for (i = 0; i < m + 1; i++) row[i] = aa->i[i] + 1;
151: for (i = 0; i < nz; i++) col[i] = aa->j[i] + 1;
152: *r = (INT_TYPE *)row;
153: *c = (INT_TYPE *)col;
154: *nnz = (INT_TYPE)nz;
155: }
156: *free = PETSC_TRUE;
157: } else {
158: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "This should not happen");
159: }
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode MatMKLPardiso_Convert_seqaij(Mat A, PetscBool sym, MatReuse reuse, PetscBool *free, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, PetscScalar **v)
164: {
165: Mat_SeqAIJ *aa = (Mat_SeqAIJ *)A->data;
166: PetscScalar *aav;
168: PetscFunctionBegin;
169: PetscCall(MatSeqAIJGetArrayRead(A, (const PetscScalar **)&aav));
170: if (!sym) { /* already in the correct format */
171: *v = aav;
172: *r = (INT_TYPE *)aa->i;
173: *c = (INT_TYPE *)aa->j;
174: *nnz = (INT_TYPE)aa->nz;
175: *free = PETSC_FALSE;
176: } else if (reuse == MAT_INITIAL_MATRIX) { /* need to get the triangular part */
177: PetscScalar *vals, *vv;
178: PetscInt *row, *col, *jj;
179: PetscInt m = A->rmap->n, nz, i;
181: nz = 0;
182: for (i = 0; i < m; i++) nz += aa->i[i + 1] - aa->diag[i];
183: PetscCall(PetscMalloc2(m + 1, &row, nz, &col));
184: PetscCall(PetscMalloc1(nz, &vals));
185: jj = col;
186: vv = vals;
188: row[0] = 0;
189: for (i = 0; i < m; i++) {
190: PetscInt *aj = aa->j + aa->diag[i];
191: PetscScalar *av = aav + aa->diag[i];
192: PetscInt rl = aa->i[i + 1] - aa->diag[i], j;
194: for (j = 0; j < rl; j++) {
195: *jj = *aj;
196: jj++;
197: aj++;
198: *vv = *av;
199: vv++;
200: av++;
201: }
202: row[i + 1] = row[i] + rl;
203: }
204: *v = vals;
205: *r = (INT_TYPE *)row;
206: *c = (INT_TYPE *)col;
207: *nnz = (INT_TYPE)nz;
208: *free = PETSC_TRUE;
209: } else {
210: PetscScalar *vv;
211: PetscInt m = A->rmap->n, i;
213: vv = *v;
214: for (i = 0; i < m; i++) {
215: PetscScalar *av = aav + aa->diag[i];
216: PetscInt rl = aa->i[i + 1] - aa->diag[i], j;
217: for (j = 0; j < rl; j++) {
218: *vv = *av;
219: vv++;
220: av++;
221: }
222: }
223: *free = PETSC_TRUE;
224: }
225: PetscCall(MatSeqAIJRestoreArrayRead(A, (const PetscScalar **)&aav));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode MatMKLPardisoSolveSchur_Private(Mat F, PetscScalar *B, PetscScalar *X)
230: {
231: Mat_MKL_PARDISO *mpardiso = (Mat_MKL_PARDISO *)F->data;
232: Mat S, Xmat, Bmat;
233: MatFactorSchurStatus schurstatus;
235: PetscFunctionBegin;
236: PetscCall(MatFactorGetSchurComplement(F, &S, &schurstatus));
237: PetscCheck(X != B || schurstatus != MAT_FACTOR_SCHUR_INVERTED, PETSC_COMM_SELF, PETSC_ERR_SUP, "X and B cannot point to the same address");
238: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, mpardiso->schur_size, mpardiso->nrhs, B, &Bmat));
239: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, mpardiso->schur_size, mpardiso->nrhs, X, &Xmat));
240: PetscCall(MatSetType(Bmat, ((PetscObject)S)->type_name));
241: PetscCall(MatSetType(Xmat, ((PetscObject)S)->type_name));
242: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
243: PetscCall(MatBindToCPU(Xmat, S->boundtocpu));
244: PetscCall(MatBindToCPU(Bmat, S->boundtocpu));
245: #endif
247: #if defined(PETSC_USE_COMPLEX)
248: PetscCheck(mpardiso->iparm[12 - 1] != 1, PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Hermitian solve not implemented yet");
249: #endif
251: switch (schurstatus) {
252: case MAT_FACTOR_SCHUR_FACTORED:
253: if (!mpardiso->iparm[12 - 1]) {
254: PetscCall(MatMatSolve(S, Bmat, Xmat));
255: } else { /* transpose solve */
256: PetscCall(MatMatSolveTranspose(S, Bmat, Xmat));
257: }
258: break;
259: case MAT_FACTOR_SCHUR_INVERTED:
260: PetscCall(MatProductCreateWithMat(S, Bmat, NULL, Xmat));
261: if (!mpardiso->iparm[12 - 1]) {
262: PetscCall(MatProductSetType(Xmat, MATPRODUCT_AB));
263: } else { /* transpose solve */
264: PetscCall(MatProductSetType(Xmat, MATPRODUCT_AtB));
265: }
266: PetscCall(MatProductSetFromOptions(Xmat));
267: PetscCall(MatProductSymbolic(Xmat));
268: PetscCall(MatProductNumeric(Xmat));
269: PetscCall(MatProductClear(Xmat));
270: break;
271: default:
272: SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Unhandled MatFactorSchurStatus %" PetscInt_FMT, F->schur_status);
273: break;
274: }
275: PetscCall(MatFactorRestoreSchurComplement(F, &S, schurstatus));
276: PetscCall(MatDestroy(&Bmat));
277: PetscCall(MatDestroy(&Xmat));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode MatFactorSetSchurIS_MKL_PARDISO(Mat F, IS is)
282: {
283: Mat_MKL_PARDISO *mpardiso = (Mat_MKL_PARDISO *)F->data;
284: const PetscScalar *arr;
285: const PetscInt *idxs;
286: PetscInt size, i;
287: PetscMPIInt csize;
288: PetscBool sorted;
290: PetscFunctionBegin;
291: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)F), &csize));
292: PetscCheck(csize <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "MKL PARDISO parallel Schur complements not yet supported from PETSc");
293: PetscCall(ISSorted(is, &sorted));
294: PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS for MKL PARDISO Schur complements needs to be sorted");
295: PetscCall(ISGetLocalSize(is, &size));
296: PetscCall(PetscFree(mpardiso->schur_work));
297: PetscCall(PetscBLASIntCast(PetscMax(mpardiso->n, 2 * size), &mpardiso->schur_work_size));
298: PetscCall(PetscMalloc1(mpardiso->schur_work_size, &mpardiso->schur_work));
299: PetscCall(MatDestroy(&F->schur));
300: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, size, size, NULL, &F->schur));
301: PetscCall(MatDenseGetArrayRead(F->schur, &arr));
302: mpardiso->schur = (PetscScalar *)arr;
303: mpardiso->schur_size = size;
304: PetscCall(MatDenseRestoreArrayRead(F->schur, &arr));
305: if (mpardiso->mtype == 2) PetscCall(MatSetOption(F->schur, MAT_SPD, PETSC_TRUE));
307: PetscCall(PetscFree(mpardiso->schur_idxs));
308: PetscCall(PetscMalloc1(size, &mpardiso->schur_idxs));
309: PetscCall(PetscArrayzero(mpardiso->perm, mpardiso->n));
310: PetscCall(ISGetIndices(is, &idxs));
311: PetscCall(PetscArraycpy(mpardiso->schur_idxs, idxs, size));
312: for (i = 0; i < size; i++) mpardiso->perm[idxs[i]] = 1;
313: PetscCall(ISRestoreIndices(is, &idxs));
314: if (size) { /* turn on Schur switch if the set of indices is not empty */
315: mpardiso->iparm[36 - 1] = 2;
316: }
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: static PetscErrorCode MatDestroy_MKL_PARDISO(Mat A)
321: {
322: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)A->data;
324: PetscFunctionBegin;
325: if (mat_mkl_pardiso->CleanUp) {
326: mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
328: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, NULL, NULL, NULL, NULL, &mat_mkl_pardiso->nrhs, mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, NULL, NULL,
329: &mat_mkl_pardiso->err);
330: }
331: PetscCall(PetscFree(mat_mkl_pardiso->perm));
332: PetscCall(PetscFree(mat_mkl_pardiso->schur_work));
333: PetscCall(PetscFree(mat_mkl_pardiso->schur_idxs));
334: if (mat_mkl_pardiso->freeaij) {
335: PetscCall(PetscFree2(mat_mkl_pardiso->ia, mat_mkl_pardiso->ja));
336: if (mat_mkl_pardiso->iparm[34] == 1) PetscCall(PetscFree(mat_mkl_pardiso->a));
337: }
338: PetscCall(PetscFree(A->data));
340: /* clear composed functions */
341: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
342: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorSetSchurIS_C", NULL));
343: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMkl_PardisoSetCntl_C", NULL));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: static PetscErrorCode MatMKLPardisoScatterSchur_Private(Mat_MKL_PARDISO *mpardiso, PetscScalar *whole, PetscScalar *schur, PetscBool reduce)
348: {
349: PetscFunctionBegin;
350: if (reduce) { /* data given for the whole matrix */
351: PetscInt i, m = 0, p = 0;
352: for (i = 0; i < mpardiso->nrhs; i++) {
353: PetscInt j;
354: for (j = 0; j < mpardiso->schur_size; j++) schur[p + j] = whole[m + mpardiso->schur_idxs[j]];
355: m += mpardiso->n;
356: p += mpardiso->schur_size;
357: }
358: } else { /* from Schur to whole */
359: PetscInt i, m = 0, p = 0;
360: for (i = 0; i < mpardiso->nrhs; i++) {
361: PetscInt j;
362: for (j = 0; j < mpardiso->schur_size; j++) whole[m + mpardiso->schur_idxs[j]] = schur[p + j];
363: m += mpardiso->n;
364: p += mpardiso->schur_size;
365: }
366: }
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode MatSolve_MKL_PARDISO(Mat A, Vec b, Vec x)
371: {
372: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)A->data;
373: PetscScalar *xarray;
374: const PetscScalar *barray;
376: PetscFunctionBegin;
377: mat_mkl_pardiso->nrhs = 1;
378: PetscCall(VecGetArrayWrite(x, &xarray));
379: PetscCall(VecGetArrayRead(b, &barray));
381: if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
382: else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
384: if (barray == xarray) { /* if the two vectors share the same memory */
385: PetscScalar *work;
386: if (!mat_mkl_pardiso->schur_work) {
387: PetscCall(PetscMalloc1(mat_mkl_pardiso->n, &work));
388: } else {
389: work = mat_mkl_pardiso->schur_work;
390: }
391: mat_mkl_pardiso->iparm[6 - 1] = 1;
392: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, mat_mkl_pardiso->a, mat_mkl_pardiso->ia, mat_mkl_pardiso->ja, NULL, &mat_mkl_pardiso->nrhs,
393: mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, (void *)xarray, (void *)work, &mat_mkl_pardiso->err);
394: if (!mat_mkl_pardiso->schur_work) PetscCall(PetscFree(work));
395: } else {
396: mat_mkl_pardiso->iparm[6 - 1] = 0;
397: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, mat_mkl_pardiso->a, mat_mkl_pardiso->ia, mat_mkl_pardiso->ja, mat_mkl_pardiso->perm,
398: &mat_mkl_pardiso->nrhs, mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_pardiso->err);
399: }
400: PetscCall(VecRestoreArrayRead(b, &barray));
402: PetscCheck(mat_mkl_pardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL PARDISO: err=%d. Please check manual", mat_mkl_pardiso->err);
404: if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
405: if (!mat_mkl_pardiso->solve_interior) {
406: PetscInt shift = mat_mkl_pardiso->schur_size;
408: PetscCall(MatFactorFactorizeSchurComplement(A));
409: /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
410: if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;
412: /* solve Schur complement */
413: PetscCall(MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso, xarray, mat_mkl_pardiso->schur_work, PETSC_TRUE));
414: PetscCall(MatMKLPardisoSolveSchur_Private(A, mat_mkl_pardiso->schur_work, mat_mkl_pardiso->schur_work + shift));
415: PetscCall(MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso, xarray, mat_mkl_pardiso->schur_work + shift, PETSC_FALSE));
416: } else { /* if we are solving for the interior problem, any value in barray[schur] forward-substituted to xarray[schur] will be neglected */
417: PetscInt i;
418: for (i = 0; i < mat_mkl_pardiso->schur_size; i++) xarray[mat_mkl_pardiso->schur_idxs[i]] = 0.;
419: }
421: /* expansion phase */
422: mat_mkl_pardiso->iparm[6 - 1] = 1;
423: mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
424: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, mat_mkl_pardiso->a, mat_mkl_pardiso->ia, mat_mkl_pardiso->ja, mat_mkl_pardiso->perm,
425: &mat_mkl_pardiso->nrhs, mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, (void *)xarray, (void *)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
426: &mat_mkl_pardiso->err);
428: PetscCheck(mat_mkl_pardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL PARDISO: err=%d. Please check manual", mat_mkl_pardiso->err);
429: mat_mkl_pardiso->iparm[6 - 1] = 0;
430: }
431: PetscCall(VecRestoreArrayWrite(x, &xarray));
432: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: static PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A, Vec b, Vec x)
437: {
438: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)A->data;
439: PetscInt oiparm12;
441: PetscFunctionBegin;
442: oiparm12 = mat_mkl_pardiso->iparm[12 - 1];
443: mat_mkl_pardiso->iparm[12 - 1] = 2;
444: PetscCall(MatSolve_MKL_PARDISO(A, b, x));
445: mat_mkl_pardiso->iparm[12 - 1] = oiparm12;
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A, Mat B, Mat X)
450: {
451: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)(A)->data;
452: const PetscScalar *barray;
453: PetscScalar *xarray;
454: PetscBool flg;
456: PetscFunctionBegin;
457: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQDENSE, &flg));
458: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix B must be MATSEQDENSE matrix");
459: if (X != B) {
460: PetscCall(PetscObjectBaseTypeCompare((PetscObject)X, MATSEQDENSE, &flg));
461: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix X must be MATSEQDENSE matrix");
462: }
464: PetscCall(MatGetSize(B, NULL, (PetscInt *)&mat_mkl_pardiso->nrhs));
466: if (mat_mkl_pardiso->nrhs > 0) {
467: PetscCall(MatDenseGetArrayRead(B, &barray));
468: PetscCall(MatDenseGetArrayWrite(X, &xarray));
470: PetscCheck(barray != xarray, PETSC_COMM_SELF, PETSC_ERR_SUP, "B and X cannot share the same memory location");
471: if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
472: else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
474: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, mat_mkl_pardiso->a, mat_mkl_pardiso->ia, mat_mkl_pardiso->ja, mat_mkl_pardiso->perm,
475: &mat_mkl_pardiso->nrhs, mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_pardiso->err);
476: PetscCheck(mat_mkl_pardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL PARDISO: err=%d. Please check manual", mat_mkl_pardiso->err);
478: PetscCall(MatDenseRestoreArrayRead(B, &barray));
479: if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */
480: PetscScalar *o_schur_work = NULL;
482: /* solve Schur complement */
483: if (!mat_mkl_pardiso->solve_interior) {
484: PetscInt shift = mat_mkl_pardiso->schur_size * mat_mkl_pardiso->nrhs, scale;
485: PetscInt mem = mat_mkl_pardiso->n * mat_mkl_pardiso->nrhs;
487: PetscCall(MatFactorFactorizeSchurComplement(A));
488: /* allocate extra memory if it is needed */
489: scale = 1;
490: if (A->schur_status == MAT_FACTOR_SCHUR_INVERTED) scale = 2;
491: mem *= scale;
492: if (mem > mat_mkl_pardiso->schur_work_size) {
493: o_schur_work = mat_mkl_pardiso->schur_work;
494: PetscCall(PetscMalloc1(mem, &mat_mkl_pardiso->schur_work));
495: }
496: /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */
497: if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0;
498: PetscCall(MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso, xarray, mat_mkl_pardiso->schur_work, PETSC_TRUE));
499: PetscCall(MatMKLPardisoSolveSchur_Private(A, mat_mkl_pardiso->schur_work, mat_mkl_pardiso->schur_work + shift));
500: PetscCall(MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso, xarray, mat_mkl_pardiso->schur_work + shift, PETSC_FALSE));
501: } else { /* if we are solving for the interior problem, any value in barray[schur,n] forward-substituted to xarray[schur,n] will be neglected */
502: PetscInt i, n, m = 0;
503: for (n = 0; n < mat_mkl_pardiso->nrhs; n++) {
504: for (i = 0; i < mat_mkl_pardiso->schur_size; i++) xarray[mat_mkl_pardiso->schur_idxs[i] + m] = 0.;
505: m += mat_mkl_pardiso->n;
506: }
507: }
509: /* expansion phase */
510: mat_mkl_pardiso->iparm[6 - 1] = 1;
511: mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
512: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, mat_mkl_pardiso->a, mat_mkl_pardiso->ia, mat_mkl_pardiso->ja, mat_mkl_pardiso->perm,
513: &mat_mkl_pardiso->nrhs, mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, (void *)xarray, (void *)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */
514: &mat_mkl_pardiso->err);
515: if (o_schur_work) { /* restore original Schur_work (minimal size) */
516: PetscCall(PetscFree(mat_mkl_pardiso->schur_work));
517: mat_mkl_pardiso->schur_work = o_schur_work;
518: }
519: PetscCheck(mat_mkl_pardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL PARDISO: err=%d. Please check manual", mat_mkl_pardiso->err);
520: mat_mkl_pardiso->iparm[6 - 1] = 0;
521: }
522: PetscCall(MatDenseRestoreArrayWrite(X, &xarray));
523: }
524: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: static PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F, Mat A, const MatFactorInfo *info)
529: {
530: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)(F)->data;
532: PetscFunctionBegin;
533: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
534: PetscCall((*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));
536: mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION;
537: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, mat_mkl_pardiso->a, mat_mkl_pardiso->ia, mat_mkl_pardiso->ja, mat_mkl_pardiso->perm,
538: &mat_mkl_pardiso->nrhs, mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, NULL, (void *)mat_mkl_pardiso->schur, &mat_mkl_pardiso->err);
539: PetscCheck(mat_mkl_pardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL PARDISO: err=%d. Please check manual", mat_mkl_pardiso->err);
541: /* report flops */
542: if (mat_mkl_pardiso->iparm[18] > 0) PetscCall(PetscLogFlops(PetscPowRealInt(10., 6) * mat_mkl_pardiso->iparm[18]));
544: if (F->schur) { /* schur output from pardiso is in row major format */
545: #if defined(PETSC_HAVE_CUDA)
546: F->schur->offloadmask = PETSC_OFFLOAD_CPU;
547: #endif
548: PetscCall(MatFactorRestoreSchurComplement(F, NULL, MAT_FACTOR_SCHUR_UNFACTORED));
549: PetscCall(MatTranspose(F->schur, MAT_INPLACE_MATRIX, &F->schur));
550: }
551: mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN;
552: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static PetscErrorCode MatSetFromOptions_MKL_PARDISO(Mat F, Mat A)
557: {
558: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)F->data;
559: PetscInt icntl, bs, threads = 1;
560: PetscBool flg;
562: PetscFunctionBegin;
563: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "MKL_PARDISO Options", "Mat");
565: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_65", "Suggested number of threads to use within MKL PARDISO", "None", threads, &threads, &flg));
566: if (flg) PetscSetMKL_PARDISOThreads((int)threads);
568: PetscCall(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));
569: if (flg) mat_mkl_pardiso->maxfct = icntl;
571: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_67", "Indicates the actual matrix for the solution phase", "None", mat_mkl_pardiso->mnum, &icntl, &flg));
572: if (flg) mat_mkl_pardiso->mnum = icntl;
574: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_68", "Message level information", "None", mat_mkl_pardiso->msglvl, &icntl, &flg));
575: if (flg) mat_mkl_pardiso->msglvl = icntl;
577: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_69", "Defines the matrix type", "None", mat_mkl_pardiso->mtype, &icntl, &flg));
578: if (flg) {
579: void *pt[IPARM_SIZE];
580: mat_mkl_pardiso->mtype = icntl;
581: icntl = mat_mkl_pardiso->iparm[34];
582: bs = mat_mkl_pardiso->iparm[36];
583: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
584: #if defined(PETSC_USE_REAL_SINGLE)
585: mat_mkl_pardiso->iparm[27] = 1;
586: #else
587: mat_mkl_pardiso->iparm[27] = 0;
588: #endif
589: mat_mkl_pardiso->iparm[34] = icntl;
590: mat_mkl_pardiso->iparm[36] = bs;
591: }
593: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_1", "Use default values (if 0)", "None", mat_mkl_pardiso->iparm[0], &icntl, &flg));
594: if (flg) mat_mkl_pardiso->iparm[0] = icntl;
596: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_2", "Fill-in reducing ordering for the input matrix", "None", mat_mkl_pardiso->iparm[1], &icntl, &flg));
597: if (flg) mat_mkl_pardiso->iparm[1] = icntl;
599: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_4", "Preconditioned CGS/CG", "None", mat_mkl_pardiso->iparm[3], &icntl, &flg));
600: if (flg) mat_mkl_pardiso->iparm[3] = icntl;
602: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_5", "User permutation", "None", mat_mkl_pardiso->iparm[4], &icntl, &flg));
603: if (flg) mat_mkl_pardiso->iparm[4] = icntl;
605: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_6", "Write solution on x", "None", mat_mkl_pardiso->iparm[5], &icntl, &flg));
606: if (flg) mat_mkl_pardiso->iparm[5] = icntl;
608: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_8", "Iterative refinement step", "None", mat_mkl_pardiso->iparm[7], &icntl, &flg));
609: if (flg) mat_mkl_pardiso->iparm[7] = icntl;
611: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_10", "Pivoting perturbation", "None", mat_mkl_pardiso->iparm[9], &icntl, &flg));
612: if (flg) mat_mkl_pardiso->iparm[9] = icntl;
614: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_11", "Scaling vectors", "None", mat_mkl_pardiso->iparm[10], &icntl, &flg));
615: if (flg) mat_mkl_pardiso->iparm[10] = icntl;
617: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_12", "Solve with transposed or conjugate transposed matrix A", "None", mat_mkl_pardiso->iparm[11], &icntl, &flg));
618: if (flg) mat_mkl_pardiso->iparm[11] = icntl;
620: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_13", "Improved accuracy using (non-) symmetric weighted matching", "None", mat_mkl_pardiso->iparm[12], &icntl, &flg));
621: if (flg) mat_mkl_pardiso->iparm[12] = icntl;
623: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_18", "Numbers of non-zero elements", "None", mat_mkl_pardiso->iparm[17], &icntl, &flg));
624: if (flg) mat_mkl_pardiso->iparm[17] = icntl;
626: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_19", "Report number of floating point operations (0 to disable)", "None", mat_mkl_pardiso->iparm[18], &icntl, &flg));
627: if (flg) mat_mkl_pardiso->iparm[18] = icntl;
629: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_21", "Pivoting for symmetric indefinite matrices", "None", mat_mkl_pardiso->iparm[20], &icntl, &flg));
630: if (flg) mat_mkl_pardiso->iparm[20] = icntl;
632: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_24", "Parallel factorization control", "None", mat_mkl_pardiso->iparm[23], &icntl, &flg));
633: if (flg) mat_mkl_pardiso->iparm[23] = icntl;
635: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_25", "Parallel forward/backward solve control", "None", mat_mkl_pardiso->iparm[24], &icntl, &flg));
636: if (flg) mat_mkl_pardiso->iparm[24] = icntl;
638: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_27", "Matrix checker", "None", mat_mkl_pardiso->iparm[26], &icntl, &flg));
639: if (flg) mat_mkl_pardiso->iparm[26] = icntl;
641: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_31", "Partial solve and computing selected components of the solution vectors", "None", mat_mkl_pardiso->iparm[30], &icntl, &flg));
642: if (flg) mat_mkl_pardiso->iparm[30] = icntl;
644: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_34", "Optimal number of threads for conditional numerical reproducibility (CNR) mode", "None", mat_mkl_pardiso->iparm[33], &icntl, &flg));
645: if (flg) mat_mkl_pardiso->iparm[33] = icntl;
647: PetscCall(PetscOptionsInt("-mat_mkl_pardiso_60", "Intel MKL PARDISO mode", "None", mat_mkl_pardiso->iparm[59], &icntl, &flg));
648: if (flg) mat_mkl_pardiso->iparm[59] = icntl;
649: PetscOptionsEnd();
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: static PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso)
654: {
655: PetscInt i, bs;
656: PetscBool match;
658: PetscFunctionBegin;
659: for (i = 0; i < IPARM_SIZE; i++) mat_mkl_pardiso->iparm[i] = 0;
660: for (i = 0; i < IPARM_SIZE; i++) mat_mkl_pardiso->pt[i] = 0;
661: #if defined(PETSC_USE_REAL_SINGLE)
662: mat_mkl_pardiso->iparm[27] = 1;
663: #else
664: mat_mkl_pardiso->iparm[27] = 0;
665: #endif
666: /* Default options for both sym and unsym */
667: mat_mkl_pardiso->iparm[0] = 1; /* Solver default parameters overridden with provided by iparm */
668: mat_mkl_pardiso->iparm[1] = 2; /* Metis reordering */
669: mat_mkl_pardiso->iparm[5] = 0; /* Write solution into x */
670: mat_mkl_pardiso->iparm[7] = 0; /* Max number of iterative refinement steps */
671: mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
672: mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
673: #if 0
674: mat_mkl_pardiso->iparm[23] = 1; /* Parallel factorization control*/
675: #endif
676: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &match, MATSEQBAIJ, MATSEQSBAIJ, ""));
677: PetscCall(MatGetBlockSize(A, &bs));
678: if (!match || bs == 1) {
679: mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */
680: mat_mkl_pardiso->n = A->rmap->N;
681: } else {
682: mat_mkl_pardiso->iparm[34] = 0; /* Cluster Sparse Solver use Fortran-style indexing for ia and ja arrays */
683: mat_mkl_pardiso->iparm[36] = bs;
684: mat_mkl_pardiso->n = A->rmap->N / bs;
685: }
686: mat_mkl_pardiso->iparm[39] = 0; /* Input: matrix/rhs/solution stored on rank-0 */
688: mat_mkl_pardiso->CleanUp = PETSC_FALSE;
689: mat_mkl_pardiso->maxfct = 1; /* Maximum number of numerical factorizations. */
690: mat_mkl_pardiso->mnum = 1; /* Which factorization to use. */
691: mat_mkl_pardiso->msglvl = 0; /* 0: do not print 1: Print statistical information in file */
692: mat_mkl_pardiso->phase = -1;
693: mat_mkl_pardiso->err = 0;
695: mat_mkl_pardiso->nrhs = 1;
696: mat_mkl_pardiso->err = 0;
697: mat_mkl_pardiso->phase = -1;
699: if (ftype == MAT_FACTOR_LU) {
700: mat_mkl_pardiso->iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
701: mat_mkl_pardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
702: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
703: } else {
704: mat_mkl_pardiso->iparm[9] = 8; /* Perturb the pivot elements with 1E-8 */
705: mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */
706: mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
707: #if defined(PETSC_USE_DEBUG)
708: mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */
709: #endif
710: }
711: PetscCall(PetscCalloc1(A->rmap->N * sizeof(INT_TYPE), &mat_mkl_pardiso->perm));
712: mat_mkl_pardiso->schur_size = 0;
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: static PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F, Mat A, const MatFactorInfo *info)
717: {
718: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)F->data;
720: PetscFunctionBegin;
721: mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
722: PetscCall(MatSetFromOptions_MKL_PARDISO(F, A));
723: /* throw away any previously computed structure */
724: if (mat_mkl_pardiso->freeaij) {
725: PetscCall(PetscFree2(mat_mkl_pardiso->ia, mat_mkl_pardiso->ja));
726: if (mat_mkl_pardiso->iparm[34] == 1) PetscCall(PetscFree(mat_mkl_pardiso->a));
727: }
728: PetscCall((*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));
729: if (mat_mkl_pardiso->iparm[34] == 1) mat_mkl_pardiso->n = A->rmap->N;
730: else mat_mkl_pardiso->n = A->rmap->N / A->rmap->bs;
732: mat_mkl_pardiso->phase = JOB_ANALYSIS;
734: /* reset flops counting if requested */
735: if (mat_mkl_pardiso->iparm[18]) mat_mkl_pardiso->iparm[18] = -1;
737: MKL_PARDISO(mat_mkl_pardiso->pt, &mat_mkl_pardiso->maxfct, &mat_mkl_pardiso->mnum, &mat_mkl_pardiso->mtype, &mat_mkl_pardiso->phase, &mat_mkl_pardiso->n, mat_mkl_pardiso->a, mat_mkl_pardiso->ia, mat_mkl_pardiso->ja, mat_mkl_pardiso->perm,
738: &mat_mkl_pardiso->nrhs, mat_mkl_pardiso->iparm, &mat_mkl_pardiso->msglvl, NULL, NULL, &mat_mkl_pardiso->err);
739: PetscCheck(mat_mkl_pardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL PARDISO: err=%d. Please check manual", mat_mkl_pardiso->err);
741: mat_mkl_pardiso->CleanUp = PETSC_TRUE;
743: if (F->factortype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO;
744: else F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO;
746: F->ops->solve = MatSolve_MKL_PARDISO;
747: F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO;
748: F->ops->matsolve = MatMatSolve_MKL_PARDISO;
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: static PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
753: {
754: PetscFunctionBegin;
755: PetscCall(MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: #if !defined(PETSC_USE_COMPLEX)
760: static PetscErrorCode MatGetInertia_MKL_PARDISO(Mat F, PetscInt *nneg, PetscInt *nzero, PetscInt *npos)
761: {
762: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)F->data;
764: PetscFunctionBegin;
765: if (nneg) *nneg = mat_mkl_pardiso->iparm[22];
766: if (npos) *npos = mat_mkl_pardiso->iparm[21];
767: if (nzero) *nzero = F->rmap->N - (mat_mkl_pardiso->iparm[22] + mat_mkl_pardiso->iparm[21]);
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
770: #endif
772: static PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F, Mat A, IS r, const MatFactorInfo *info)
773: {
774: PetscFunctionBegin;
775: PetscCall(MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info));
776: F->ops->getinertia = NULL;
777: #if !defined(PETSC_USE_COMPLEX)
778: F->ops->getinertia = MatGetInertia_MKL_PARDISO;
779: #endif
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: static PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer)
784: {
785: PetscBool iascii;
786: PetscViewerFormat format;
787: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)A->data;
788: PetscInt i;
790: PetscFunctionBegin;
791: if (A->ops->solve != MatSolve_MKL_PARDISO) PetscFunctionReturn(PETSC_SUCCESS);
793: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
794: if (iascii) {
795: PetscCall(PetscViewerGetFormat(viewer, &format));
796: if (format == PETSC_VIEWER_ASCII_INFO) {
797: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO run parameters:\n"));
798: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO phase: %d \n", mat_mkl_pardiso->phase));
799: for (i = 1; i <= 64; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO iparm[%d]: %d \n", i, mat_mkl_pardiso->iparm[i - 1]));
800: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct));
801: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO mnum: %d \n", mat_mkl_pardiso->mnum));
802: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO mtype: %d \n", mat_mkl_pardiso->mtype));
803: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO n: %d \n", mat_mkl_pardiso->n));
804: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs));
805: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl));
806: }
807: }
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info)
812: {
813: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)A->data;
815: PetscFunctionBegin;
816: info->block_size = 1.0;
817: info->nz_used = mat_mkl_pardiso->iparm[17];
818: info->nz_allocated = mat_mkl_pardiso->iparm[17];
819: info->nz_unneeded = 0.0;
820: info->assemblies = 0.0;
821: info->mallocs = 0.0;
822: info->memory = 0.0;
823: info->fill_ratio_given = 0;
824: info->fill_ratio_needed = 0;
825: info->factor_mallocs = 0;
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }
829: static PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F, PetscInt icntl, PetscInt ival)
830: {
831: PetscInt backup, bs;
832: Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO *)F->data;
834: PetscFunctionBegin;
835: if (icntl <= 64) {
836: mat_mkl_pardiso->iparm[icntl - 1] = ival;
837: } else {
838: if (icntl == 65) PetscSetMKL_PARDISOThreads(ival);
839: else if (icntl == 66) mat_mkl_pardiso->maxfct = ival;
840: else if (icntl == 67) mat_mkl_pardiso->mnum = ival;
841: else if (icntl == 68) mat_mkl_pardiso->msglvl = ival;
842: else if (icntl == 69) {
843: void *pt[IPARM_SIZE];
844: backup = mat_mkl_pardiso->iparm[34];
845: bs = mat_mkl_pardiso->iparm[36];
846: mat_mkl_pardiso->mtype = ival;
847: MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm);
848: #if defined(PETSC_USE_REAL_SINGLE)
849: mat_mkl_pardiso->iparm[27] = 1;
850: #else
851: mat_mkl_pardiso->iparm[27] = 0;
852: #endif
853: mat_mkl_pardiso->iparm[34] = backup;
854: mat_mkl_pardiso->iparm[36] = bs;
855: } else if (icntl == 70) mat_mkl_pardiso->solve_interior = (PetscBool) !!ival;
856: }
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*@
861: MatMkl_PardisoSetCntl - Set MKL PARDISO <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html> parameters
863: Logically Collective
865: Input Parameters:
866: + F - the factored matrix obtained by calling `MatGetFactor()`
867: . icntl - index of MKL PARDISO parameter
868: - ival - value of MKL PARDISO parameter
870: Options Database Key:
871: . -mat_mkl_pardiso_<icntl> <ival> - change the option numbered icntl to the value ival
873: Level: beginner
875: .seealso: [](ch_matrices), `Mat`, `MATSOLVERMKL_PARDISO`, `MatGetFactor()`
876: @*/
877: PetscErrorCode MatMkl_PardisoSetCntl(Mat F, PetscInt icntl, PetscInt ival)
878: {
879: PetscFunctionBegin;
880: PetscTryMethod(F, "MatMkl_PardisoSetCntl_C", (Mat, PetscInt, PetscInt), (F, icntl, ival));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: /*MC
885: MATSOLVERMKL_PARDISO - A matrix type providing direct solvers, LU, for
886: `MATSEQAIJ` matrices via the external package MKL PARDISO
887: <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/onemkl-pardiso-parallel-direct-sparse-solver-iface.html>.
889: Use `-pc_type lu` `-pc_factor_mat_solver_type mkl_pardiso` to use this direct solver
891: Options Database Keys:
892: + -mat_mkl_pardiso_65 - Suggested number of threads to use within MKL PARDISO
893: . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
894: . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase
895: . -mat_mkl_pardiso_68 - Message level information, use 1 to get detailed information on the solver options
896: . -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
897: . -mat_mkl_pardiso_1 - Use default values
898: . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix
899: . -mat_mkl_pardiso_4 - Preconditioned CGS/CG
900: . -mat_mkl_pardiso_5 - User permutation
901: . -mat_mkl_pardiso_6 - Write solution on x
902: . -mat_mkl_pardiso_8 - Iterative refinement step
903: . -mat_mkl_pardiso_10 - Pivoting perturbation
904: . -mat_mkl_pardiso_11 - Scaling vectors
905: . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A
906: . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching
907: . -mat_mkl_pardiso_18 - Numbers of non-zero elements
908: . -mat_mkl_pardiso_19 - Report number of floating point operations
909: . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices
910: . -mat_mkl_pardiso_24 - Parallel factorization control
911: . -mat_mkl_pardiso_25 - Parallel forward/backward solve control
912: . -mat_mkl_pardiso_27 - Matrix checker
913: . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors
914: . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
915: - -mat_mkl_pardiso_60 - Intel MKL PARDISO mode
917: Level: beginner
919: Notes:
920: Use `-mat_mkl_pardiso_68 1` to display the number of threads the solver is using. MKL does not provide a way to directly access this
921: information.
923: For more information on the options check the MKL PARDISO manual
925: .seealso: [](ch_matrices), `Mat`, `MATSEQAIJ`, `PCFactorSetMatSolverType()`, `MatSolverType`, `MatMkl_PardisoSetCntl()`, `MATSOLVERMKL_CPARDISO`
926: M*/
927: static PetscErrorCode MatFactorGetSolverType_mkl_pardiso(Mat A, MatSolverType *type)
928: {
929: PetscFunctionBegin;
930: *type = MATSOLVERMKL_PARDISO;
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A, MatFactorType ftype, Mat *F)
935: {
936: Mat B;
937: Mat_MKL_PARDISO *mat_mkl_pardiso;
938: PetscBool isSeqAIJ, isSeqBAIJ, isSeqSBAIJ;
940: PetscFunctionBegin;
941: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isSeqAIJ));
942: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQBAIJ, &isSeqBAIJ));
943: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &isSeqSBAIJ));
944: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
945: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
946: PetscCall(PetscStrallocpy("mkl_pardiso", &((PetscObject)B)->type_name));
947: PetscCall(MatSetUp(B));
949: PetscCall(PetscNew(&mat_mkl_pardiso));
950: B->data = mat_mkl_pardiso;
952: PetscCall(MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso));
953: if (ftype == MAT_FACTOR_LU) {
954: B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO;
955: B->factortype = MAT_FACTOR_LU;
956: mat_mkl_pardiso->needsym = PETSC_FALSE;
957: if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
958: else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
959: else {
960: PetscCheck(!isSeqSBAIJ, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for MKL PARDISO LU factor with SEQSBAIJ format! Use MAT_FACTOR_CHOLESKY instead");
961: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for MKL PARDISO LU with %s format", ((PetscObject)A)->type_name);
962: }
963: #if defined(PETSC_USE_COMPLEX)
964: mat_mkl_pardiso->mtype = 13;
965: #else
966: mat_mkl_pardiso->mtype = 11;
967: #endif
968: } else {
969: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO;
970: B->factortype = MAT_FACTOR_CHOLESKY;
971: if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij;
972: else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij;
973: else if (isSeqSBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqsbaij;
974: else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for PARDISO CHOLESKY with %s format", ((PetscObject)A)->type_name);
976: mat_mkl_pardiso->needsym = PETSC_TRUE;
977: #if !defined(PETSC_USE_COMPLEX)
978: if (A->spd == PETSC_BOOL3_TRUE) mat_mkl_pardiso->mtype = 2;
979: else mat_mkl_pardiso->mtype = -2;
980: #else
981: mat_mkl_pardiso->mtype = 6;
982: PetscCheck(A->hermitian != PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for MKL PARDISO CHOLESKY with Hermitian matrices! Use MAT_FACTOR_LU instead");
983: #endif
984: }
985: B->ops->destroy = MatDestroy_MKL_PARDISO;
986: B->ops->view = MatView_MKL_PARDISO;
987: B->ops->getinfo = MatGetInfo_MKL_PARDISO;
988: B->factortype = ftype;
989: B->assembled = PETSC_TRUE;
991: PetscCall(PetscFree(B->solvertype));
992: PetscCall(PetscStrallocpy(MATSOLVERMKL_PARDISO, &B->solvertype));
994: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_mkl_pardiso));
995: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorSetSchurIS_C", MatFactorSetSchurIS_MKL_PARDISO));
996: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMkl_PardisoSetCntl_C", MatMkl_PardisoSetCntl_MKL_PARDISO));
998: *F = B;
999: PetscFunctionReturn(PETSC_SUCCESS);
1000: }
1002: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void)
1003: {
1004: PetscFunctionBegin;
1005: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_PARDISO, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_mkl_pardiso));
1006: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_PARDISO, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_aij_mkl_pardiso));
1007: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_PARDISO, MATSEQBAIJ, MAT_FACTOR_LU, MatGetFactor_aij_mkl_pardiso));
1008: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_PARDISO, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_aij_mkl_pardiso));
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }