Actual source code: mkl_cpardiso.c
1: #include <../src/mat/impls/aij/mpi/mpiaij.h>
2: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
4: #if PetscDefined(HAVE_MKL_INTEL_ILP64)
5: #define MKL_ILP64
6: #endif
7: #include <mkl.h>
8: #include <mkl_cluster_sparse_solver.h>
10: /*
11: * Possible mkl_cpardiso phases that controls the execution of the solver.
12: * For more information check mkl_cpardiso manual.
13: */
14: #define JOB_ANALYSIS 11
15: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
16: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
17: #define JOB_NUMERICAL_FACTORIZATION 22
18: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
19: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
20: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
21: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
22: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
23: #define JOB_RELEASE_OF_LU_MEMORY 0
24: #define JOB_RELEASE_OF_ALL_MEMORY -1
26: #define IPARM_SIZE 64
27: #define INT_TYPE MKL_INT
29: static const char *Err_MSG_CPardiso(int errNo)
30: {
31: switch (errNo) {
32: case -1:
33: return "input inconsistent";
34: break;
35: case -2:
36: return "not enough memory";
37: break;
38: case -3:
39: return "reordering problem";
40: break;
41: case -4:
42: return "zero pivot, numerical factorization or iterative refinement problem";
43: break;
44: case -5:
45: return "unclassified (internal) error";
46: break;
47: case -6:
48: return "preordering failed (matrix types 11, 13 only)";
49: break;
50: case -7:
51: return "diagonal matrix problem";
52: break;
53: case -8:
54: return "32-bit integer overflow problem";
55: break;
56: case -9:
57: return "not enough memory for OOC";
58: break;
59: case -10:
60: return "problems with opening OOC temporary files";
61: break;
62: case -11:
63: return "read/write problems with the OOC data file";
64: break;
65: default:
66: return "unknown error";
67: }
68: }
70: #define PetscCallCluster(f) PetscCallExternalVoid("cluster_sparse_solver", f);
72: /*
73: * Internal data structure.
74: * For more information check mkl_cpardiso manual.
75: */
77: typedef struct {
78: /* Configuration vector */
79: INT_TYPE iparm[IPARM_SIZE];
81: /*
82: * Internal mkl_cpardiso memory location.
83: * After the first call to mkl_cpardiso do not modify pt, as that could cause a serious memory leak.
84: */
85: void *pt[IPARM_SIZE];
87: MPI_Fint comm_mkl_cpardiso;
89: /* Basic mkl_cpardiso info*/
90: INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;
92: /* Matrix values and matrix nonzero structure */
93: PetscScalar *a;
95: INT_TYPE *ia, *ja;
97: /* Number of non-zero elements */
98: INT_TYPE nz;
100: /* Row permutaton vector*/
101: INT_TYPE *perm;
103: /* Define is matrix preserve sparse structure. */
104: MatStructure matstruc;
106: PetscErrorCode (*ConvertToTriples)(Mat, MatReuse, PetscInt *, PetscInt **, PetscInt **, PetscScalar **);
108: /* True if mkl_cpardiso function have been used. */
109: PetscBool CleanUp;
110: } Mat_MKL_CPARDISO;
112: /*
113: * Copy the elements of matrix A.
114: * Input:
115: * - Mat A: MATSEQAIJ matrix
116: * - int shift: matrix index.
117: * - 0 for c representation
118: * - 1 for fortran representation
119: * - MatReuse reuse:
120: * - MAT_INITIAL_MATRIX: Create a new aij representation
121: * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
122: * Output:
123: * - int *nnz: Number of nonzero-elements.
124: * - int **r pointer to i index
125: * - int **c pointer to j elements
126: * - MATRIXTYPE **v: Non-zero elements
127: */
128: static PetscErrorCode MatCopy_seqaij_seqaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
129: {
130: Mat_SeqAIJ *aa = (Mat_SeqAIJ *)A->data;
132: PetscFunctionBegin;
133: *v = aa->a;
134: if (reuse == MAT_INITIAL_MATRIX) {
135: *r = (INT_TYPE *)aa->i;
136: *c = (INT_TYPE *)aa->j;
137: *nnz = aa->nz;
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: static PetscErrorCode MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
143: {
144: const PetscInt *ai, *aj, *bi, *bj, *garray, m = A->rmap->n, *ajj, *bjj;
145: PetscInt rstart, nz, i, j, countA, countB;
146: PetscInt *row, *col;
147: const PetscScalar *av, *bv;
148: PetscScalar *val;
149: Mat_MPIAIJ *mat = (Mat_MPIAIJ *)A->data;
150: Mat_SeqAIJ *aa = (Mat_SeqAIJ *)mat->A->data;
151: Mat_SeqAIJ *bb = (Mat_SeqAIJ *)mat->B->data;
152: PetscInt colA_start, jB, jcol;
154: PetscFunctionBegin;
155: ai = aa->i;
156: aj = aa->j;
157: bi = bb->i;
158: bj = bb->j;
159: rstart = A->rmap->rstart;
160: av = aa->a;
161: bv = bb->a;
163: garray = mat->garray;
165: if (reuse == MAT_INITIAL_MATRIX) {
166: nz = aa->nz + bb->nz;
167: *nnz = nz;
168: PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz, &val));
169: *r = row;
170: *c = col;
171: *v = val;
172: } else {
173: row = *r;
174: col = *c;
175: val = *v;
176: }
178: nz = 0;
179: for (i = 0; i < m; i++) {
180: row[i] = nz;
181: countA = ai[i + 1] - ai[i];
182: countB = bi[i + 1] - bi[i];
183: ajj = aj + ai[i]; /* ptr to the beginning of this row */
184: bjj = bj + bi[i];
186: /* B part, smaller col index */
187: colA_start = rstart + ajj[0]; /* the smallest global col index of A */
188: jB = 0;
189: for (j = 0; j < countB; j++) {
190: jcol = garray[bjj[j]];
191: if (jcol > colA_start) break;
192: col[nz] = jcol;
193: val[nz++] = *bv++;
194: }
195: jB = j;
197: /* A part */
198: for (j = 0; j < countA; j++) {
199: col[nz] = rstart + ajj[j];
200: val[nz++] = *av++;
201: }
203: /* B part, larger col index */
204: for (j = jB; j < countB; j++) {
205: col[nz] = garray[bjj[j]];
206: val[nz++] = *bv++;
207: }
208: }
209: row[m] = nz;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: static PetscErrorCode MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
214: {
215: const PetscInt *ai, *aj, *bi, *bj, *garray, bs = A->rmap->bs, bs2 = bs * bs, m = A->rmap->n / bs, *ajj, *bjj;
216: PetscInt rstart, nz, i, j, countA, countB;
217: PetscInt *row, *col;
218: const PetscScalar *av, *bv;
219: PetscScalar *val;
220: Mat_MPIBAIJ *mat = (Mat_MPIBAIJ *)A->data;
221: Mat_SeqBAIJ *aa = (Mat_SeqBAIJ *)mat->A->data;
222: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ *)mat->B->data;
223: PetscInt colA_start, jB, jcol;
225: PetscFunctionBegin;
226: ai = aa->i;
227: aj = aa->j;
228: bi = bb->i;
229: bj = bb->j;
230: rstart = A->rmap->rstart / bs;
231: av = aa->a;
232: bv = bb->a;
234: garray = mat->garray;
236: if (reuse == MAT_INITIAL_MATRIX) {
237: nz = aa->nz + bb->nz;
238: *nnz = nz;
239: PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz * bs2, &val));
240: *r = row;
241: *c = col;
242: *v = val;
243: } else {
244: row = *r;
245: col = *c;
246: val = *v;
247: }
249: nz = 0;
250: for (i = 0; i < m; i++) {
251: row[i] = nz + 1;
252: countA = ai[i + 1] - ai[i];
253: countB = bi[i + 1] - bi[i];
254: ajj = aj + ai[i]; /* ptr to the beginning of this row */
255: bjj = bj + bi[i];
257: /* B part, smaller col index */
258: colA_start = rstart + (countA > 0 ? ajj[0] : 0); /* the smallest global col index of A */
259: jB = 0;
260: for (j = 0; j < countB; j++) {
261: jcol = garray[bjj[j]];
262: if (jcol > colA_start) break;
263: col[nz++] = jcol + 1;
264: }
265: jB = j;
266: PetscCall(PetscArraycpy(val, bv, jB * bs2));
267: val += jB * bs2;
268: bv += jB * bs2;
270: /* A part */
271: for (j = 0; j < countA; j++) col[nz++] = rstart + ajj[j] + 1;
272: PetscCall(PetscArraycpy(val, av, countA * bs2));
273: val += countA * bs2;
274: av += countA * bs2;
276: /* B part, larger col index */
277: for (j = jB; j < countB; j++) col[nz++] = garray[bjj[j]] + 1;
278: PetscCall(PetscArraycpy(val, bv, (countB - jB) * bs2));
279: val += (countB - jB) * bs2;
280: bv += (countB - jB) * bs2;
281: }
282: row[m] = nz + 1;
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: static PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
287: {
288: const PetscInt *ai, *aj, *bi, *bj, *garray, bs = A->rmap->bs, bs2 = bs * bs, m = A->rmap->n / bs, *ajj, *bjj;
289: PetscInt rstart, nz, i, j, countA, countB;
290: PetscInt *row, *col;
291: const PetscScalar *av, *bv;
292: PetscScalar *val;
293: Mat_MPISBAIJ *mat = (Mat_MPISBAIJ *)A->data;
294: Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ *)mat->A->data;
295: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ *)mat->B->data;
297: PetscFunctionBegin;
298: ai = aa->i;
299: aj = aa->j;
300: bi = bb->i;
301: bj = bb->j;
302: rstart = A->rmap->rstart / bs;
303: av = aa->a;
304: bv = bb->a;
306: garray = mat->garray;
308: if (reuse == MAT_INITIAL_MATRIX) {
309: nz = aa->nz + bb->nz;
310: *nnz = nz;
311: PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz * bs2, &val));
312: *r = row;
313: *c = col;
314: *v = val;
315: } else {
316: row = *r;
317: col = *c;
318: val = *v;
319: }
321: nz = 0;
322: for (i = 0; i < m; i++) {
323: row[i] = nz + 1;
324: countA = ai[i + 1] - ai[i];
325: countB = bi[i + 1] - bi[i];
326: ajj = aj + ai[i]; /* ptr to the beginning of this row */
327: bjj = bj + bi[i];
329: /* A part */
330: for (j = 0; j < countA; j++) col[nz++] = rstart + ajj[j] + 1;
331: PetscCall(PetscArraycpy(val, av, countA * bs2));
332: val += countA * bs2;
333: av += countA * bs2;
335: /* B part, larger col index */
336: for (j = 0; j < countB; j++) col[nz++] = garray[bjj[j]] + 1;
337: PetscCall(PetscArraycpy(val, bv, countB * bs2));
338: val += countB * bs2;
339: bv += countB * bs2;
340: }
341: row[m] = nz + 1;
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*
346: * Free memory for Mat_MKL_CPARDISO structure and pointers to objects.
347: */
348: static PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A)
349: {
350: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
351: MPI_Comm comm;
353: PetscFunctionBegin;
354: /* Terminate instance, deallocate memories */
355: if (mat_mkl_cpardiso->CleanUp) {
356: mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
358: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, NULL, NULL, NULL, mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs,
359: mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
360: }
361: if (mat_mkl_cpardiso->ConvertToTriples != MatCopy_seqaij_seqaij_MKL_CPARDISO) PetscCall(PetscFree3(mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja, mat_mkl_cpardiso->a));
362: comm = MPI_Comm_f2c(mat_mkl_cpardiso->comm_mkl_cpardiso);
363: PetscCallMPI(MPI_Comm_free(&comm));
364: PetscCall(PetscFree(A->data));
366: /* clear composed functions */
367: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
368: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMkl_CPardisoSetCntl_C", NULL));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: /*
373: * Computes Ax = b
374: */
375: static PetscErrorCode MatSolve_MKL_CPARDISO(Mat A, Vec b, Vec x)
376: {
377: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
378: PetscScalar *xarray;
379: const PetscScalar *barray;
381: PetscFunctionBegin;
382: mat_mkl_cpardiso->nrhs = 1;
383: PetscCall(VecGetArray(x, &xarray));
384: PetscCall(VecGetArrayRead(b, &barray));
386: /* solve phase */
387: mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
388: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
389: mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
390: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
392: PetscCall(VecRestoreArray(x, &xarray));
393: PetscCall(VecRestoreArrayRead(b, &barray));
394: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: static PetscErrorCode MatForwardSolve_MKL_CPARDISO(Mat A, Vec b, Vec x)
399: {
400: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
401: PetscScalar *xarray;
402: const PetscScalar *barray;
404: PetscFunctionBegin;
405: mat_mkl_cpardiso->nrhs = 1;
406: PetscCall(VecGetArray(x, &xarray));
407: PetscCall(VecGetArrayRead(b, &barray));
409: /* solve phase */
410: mat_mkl_cpardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
411: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
412: mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
413: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
415: PetscCall(VecRestoreArray(x, &xarray));
416: PetscCall(VecRestoreArrayRead(b, &barray));
417: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: static PetscErrorCode MatBackwardSolve_MKL_CPARDISO(Mat A, Vec b, Vec x)
422: {
423: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
424: PetscScalar *xarray;
425: const PetscScalar *barray;
427: PetscFunctionBegin;
428: mat_mkl_cpardiso->nrhs = 1;
429: PetscCall(VecGetArray(x, &xarray));
430: PetscCall(VecGetArrayRead(b, &barray));
432: /* solve phase */
433: mat_mkl_cpardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
434: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
435: mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
436: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
438: PetscCall(VecRestoreArray(x, &xarray));
439: PetscCall(VecRestoreArrayRead(b, &barray));
440: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: static PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A, Vec b, Vec x)
445: {
446: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
448: PetscFunctionBegin;
449: mat_mkl_cpardiso->iparm[12 - 1] = PetscDefined(USE_COMPLEX) ? 1 : 2;
450: PetscCall(MatSolve_MKL_CPARDISO(A, b, x));
451: mat_mkl_cpardiso->iparm[12 - 1] = 0;
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: static PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A, Mat B, Mat X)
456: {
457: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
458: PetscScalar *xarray;
459: const PetscScalar *barray;
461: PetscFunctionBegin;
462: PetscCall(MatGetSize(B, NULL, (PetscInt *)&mat_mkl_cpardiso->nrhs));
464: if (mat_mkl_cpardiso->nrhs > 0) {
465: PetscCall(MatDenseGetArrayRead(B, &barray));
466: PetscCall(MatDenseGetArray(X, &xarray));
468: PetscCheck(barray != xarray, PETSC_COMM_SELF, PETSC_ERR_SUP, "B and X cannot share the same memory location");
470: /* solve phase */
471: mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
472: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
473: mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
474: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
475: PetscCall(MatDenseRestoreArrayRead(B, &barray));
476: PetscCall(MatDenseRestoreArray(X, &xarray));
477: }
478: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: /*
483: * LU Decomposition
484: */
485: static PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F, Mat A, const MatFactorInfo *info)
486: {
487: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
489: PetscFunctionBegin;
490: mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
491: PetscCall((*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_REUSE_MATRIX, &mat_mkl_cpardiso->nz, &mat_mkl_cpardiso->ia, &mat_mkl_cpardiso->ja, &mat_mkl_cpardiso->a));
493: mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION;
494: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
495: mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, &mat_mkl_cpardiso->err));
496: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
498: mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
499: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
500: PetscFunctionReturn(PETSC_SUCCESS);
501: }
503: /* Sets mkl_cpardiso options from the options database */
504: static PetscErrorCode MatSetFromOptions_MKL_CPARDISO(Mat F, Mat A)
505: {
506: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
507: PetscInt icntl, threads;
508: PetscBool flg;
510: PetscFunctionBegin;
511: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "MKL Cluster PARDISO Options", "Mat");
512: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_65", "Suggested number of threads to use within MKL Cluster PARDISO", "None", threads, &threads, &flg));
513: if (flg) mkl_set_num_threads((int)threads);
515: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_66", "Maximum number of factors with identical sparsity structure that must be kept in memory at the same time", "None", mat_mkl_cpardiso->maxfct, &icntl, &flg));
516: if (flg) mat_mkl_cpardiso->maxfct = icntl;
518: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_67", "Indicates the actual matrix for the solution phase", "None", mat_mkl_cpardiso->mnum, &icntl, &flg));
519: if (flg) mat_mkl_cpardiso->mnum = icntl;
521: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_68", "Message level information", "None", mat_mkl_cpardiso->msglvl, &icntl, &flg));
522: if (flg) mat_mkl_cpardiso->msglvl = icntl;
524: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_69", "Defines the matrix type", "None", mat_mkl_cpardiso->mtype, &icntl, &flg));
525: if (flg) mat_mkl_cpardiso->mtype = icntl;
526: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_1", "Use default values", "None", mat_mkl_cpardiso->iparm[0], &icntl, &flg));
528: if (flg && icntl != 0) {
529: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_2", "Fill-in reducing ordering for the input matrix", "None", mat_mkl_cpardiso->iparm[1], &icntl, &flg));
530: if (flg) mat_mkl_cpardiso->iparm[1] = icntl;
532: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_4", "Preconditioned CGS/CG", "None", mat_mkl_cpardiso->iparm[3], &icntl, &flg));
533: if (flg) mat_mkl_cpardiso->iparm[3] = icntl;
535: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_5", "User permutation", "None", mat_mkl_cpardiso->iparm[4], &icntl, &flg));
536: if (flg) mat_mkl_cpardiso->iparm[4] = icntl;
538: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_6", "Write solution on x", "None", mat_mkl_cpardiso->iparm[5], &icntl, &flg));
539: if (flg) mat_mkl_cpardiso->iparm[5] = icntl;
541: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_8", "Iterative refinement step", "None", mat_mkl_cpardiso->iparm[7], &icntl, &flg));
542: if (flg) mat_mkl_cpardiso->iparm[7] = icntl;
544: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_10", "Pivoting perturbation", "None", mat_mkl_cpardiso->iparm[9], &icntl, &flg));
545: if (flg) mat_mkl_cpardiso->iparm[9] = icntl;
547: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_11", "Scaling vectors", "None", mat_mkl_cpardiso->iparm[10], &icntl, &flg));
548: if (flg) mat_mkl_cpardiso->iparm[10] = icntl;
550: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_12", "Solve with transposed or conjugate transposed matrix A", "None", mat_mkl_cpardiso->iparm[11], &icntl, &flg));
551: if (flg) mat_mkl_cpardiso->iparm[11] = icntl;
553: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_13", "Improved accuracy using (non-) symmetric weighted matching", "None", mat_mkl_cpardiso->iparm[12], &icntl, &flg));
554: if (flg) mat_mkl_cpardiso->iparm[12] = icntl;
556: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_18", "Numbers of non-zero elements", "None", mat_mkl_cpardiso->iparm[17], &icntl, &flg));
557: if (flg) mat_mkl_cpardiso->iparm[17] = icntl;
559: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_19", "Report number of floating point operations", "None", mat_mkl_cpardiso->iparm[18], &icntl, &flg));
560: if (flg) mat_mkl_cpardiso->iparm[18] = icntl;
562: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_21", "Pivoting for symmetric indefinite matrices", "None", mat_mkl_cpardiso->iparm[20], &icntl, &flg));
563: if (flg) mat_mkl_cpardiso->iparm[20] = icntl;
565: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_24", "Parallel factorization control", "None", mat_mkl_cpardiso->iparm[23], &icntl, &flg));
566: if (flg) mat_mkl_cpardiso->iparm[23] = icntl;
568: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_25", "Parallel forward/backward solve control", "None", mat_mkl_cpardiso->iparm[24], &icntl, &flg));
569: if (flg) mat_mkl_cpardiso->iparm[24] = icntl;
571: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_27", "Matrix checker", "None", mat_mkl_cpardiso->iparm[26], &icntl, &flg));
572: if (flg) mat_mkl_cpardiso->iparm[26] = icntl;
574: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_31", "Partial solve and computing selected components of the solution vectors", "None", mat_mkl_cpardiso->iparm[30], &icntl, &flg));
575: if (flg) mat_mkl_cpardiso->iparm[30] = icntl;
577: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_34", "Optimal number of threads for conditional numerical reproducibility (CNR) mode", "None", mat_mkl_cpardiso->iparm[33], &icntl, &flg));
578: if (flg) mat_mkl_cpardiso->iparm[33] = icntl;
580: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_60", "Intel MKL Cluster PARDISO mode", "None", mat_mkl_cpardiso->iparm[59], &icntl, &flg));
581: if (flg) mat_mkl_cpardiso->iparm[59] = icntl;
582: }
584: PetscOptionsEnd();
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
588: static PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso)
589: {
590: PetscInt bs;
591: PetscBool match;
592: PetscMPIInt size;
593: MPI_Comm comm;
595: PetscFunctionBegin;
596: PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)A), &comm));
597: PetscCallMPI(MPI_Comm_size(comm, &size));
598: mat_mkl_cpardiso->comm_mkl_cpardiso = MPI_Comm_c2f(comm);
600: mat_mkl_cpardiso->CleanUp = PETSC_FALSE;
601: mat_mkl_cpardiso->maxfct = 1;
602: mat_mkl_cpardiso->mnum = 1;
603: mat_mkl_cpardiso->n = A->rmap->N;
604: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
605: mat_mkl_cpardiso->msglvl = 0;
606: mat_mkl_cpardiso->nrhs = 1;
607: mat_mkl_cpardiso->err = 0;
608: mat_mkl_cpardiso->phase = -1;
609: mat_mkl_cpardiso->mtype = PetscDefined(USE_COMPLEX) ? 13 : 11;
611: mat_mkl_cpardiso->iparm[27] = PetscDefined(USE_REAL_SINGLE) ? 1 : 0;
613: mat_mkl_cpardiso->iparm[0] = 1; /* Solver default parameters overridden with provided by iparm */
614: mat_mkl_cpardiso->iparm[1] = 2; /* Use METIS for fill-in reordering */
615: mat_mkl_cpardiso->iparm[5] = 0; /* Write solution into x */
616: mat_mkl_cpardiso->iparm[7] = 2; /* Max number of iterative refinement steps */
617: mat_mkl_cpardiso->iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
618: mat_mkl_cpardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
619: mat_mkl_cpardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
620: mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
621: mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
622: mat_mkl_cpardiso->iparm[26] = 1; /* Check input data for correctness */
624: mat_mkl_cpardiso->iparm[39] = 0;
625: if (size > 1) {
626: mat_mkl_cpardiso->iparm[39] = 2;
627: mat_mkl_cpardiso->iparm[40] = A->rmap->rstart;
628: mat_mkl_cpardiso->iparm[41] = A->rmap->rend - 1;
629: }
630: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &match, MATMPIBAIJ, MATMPISBAIJ, ""));
631: if (match) {
632: PetscCall(MatGetBlockSize(A, &bs));
633: mat_mkl_cpardiso->iparm[36] = bs;
634: mat_mkl_cpardiso->iparm[40] /= bs;
635: mat_mkl_cpardiso->iparm[41] /= bs;
636: mat_mkl_cpardiso->iparm[40]++;
637: mat_mkl_cpardiso->iparm[41]++;
638: mat_mkl_cpardiso->iparm[34] = 0; /* Fortran style */
639: } else {
640: mat_mkl_cpardiso->iparm[34] = 1; /* C style */
641: }
643: mat_mkl_cpardiso->perm = 0;
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: /*
648: * Symbolic decomposition. Mkl_Pardiso analysis phase.
649: */
650: static PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
651: {
652: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
654: PetscFunctionBegin;
655: mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
657: /* Set MKL_CPARDISO options from the options database */
658: PetscCall(MatSetFromOptions_MKL_CPARDISO(F, A));
659: PetscCall((*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_INITIAL_MATRIX, &mat_mkl_cpardiso->nz, &mat_mkl_cpardiso->ia, &mat_mkl_cpardiso->ja, &mat_mkl_cpardiso->a));
661: mat_mkl_cpardiso->n = A->rmap->N;
662: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
664: /* analysis phase */
665: mat_mkl_cpardiso->phase = JOB_ANALYSIS;
667: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
668: mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
669: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\".Check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
671: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
672: F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO;
673: F->ops->solve = MatSolve_MKL_CPARDISO;
674: F->ops->forwardsolve = MatForwardSolve_MKL_CPARDISO;
675: F->ops->backwardsolve = MatBackwardSolve_MKL_CPARDISO;
676: F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO;
677: F->ops->matsolve = MatMatSolve_MKL_CPARDISO;
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: static PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_CPARDISO(Mat F, Mat A, IS perm, const MatFactorInfo *info)
682: {
683: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
685: PetscFunctionBegin;
686: mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
688: /* Set MKL_CPARDISO options from the options database */
689: PetscCall(MatSetFromOptions_MKL_CPARDISO(F, A));
690: PetscCall((*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_INITIAL_MATRIX, &mat_mkl_cpardiso->nz, &mat_mkl_cpardiso->ia, &mat_mkl_cpardiso->ja, &mat_mkl_cpardiso->a));
692: mat_mkl_cpardiso->n = A->rmap->N;
693: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
694: PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead");
695: if (A->spd == PETSC_BOOL3_TRUE) mat_mkl_cpardiso->mtype = 2;
696: else mat_mkl_cpardiso->mtype = -2;
698: /* analysis phase */
699: mat_mkl_cpardiso->phase = JOB_ANALYSIS;
701: PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
702: mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
703: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\".Check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
705: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
706: F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_CPARDISO;
707: F->ops->solve = MatSolve_MKL_CPARDISO;
708: F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO;
709: F->ops->matsolve = MatMatSolve_MKL_CPARDISO;
710: if (A->spd == PETSC_BOOL3_TRUE) {
711: F->ops->forwardsolve = MatForwardSolve_MKL_CPARDISO;
712: F->ops->backwardsolve = MatBackwardSolve_MKL_CPARDISO;
713: }
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
717: static PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer)
718: {
719: PetscBool isascii;
720: PetscViewerFormat format;
721: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
723: PetscFunctionBegin;
724: /* check if matrix is mkl_cpardiso type */
725: if (A->ops->solve != MatSolve_MKL_CPARDISO) PetscFunctionReturn(PETSC_SUCCESS);
727: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
728: if (isascii) {
729: PetscCall(PetscViewerGetFormat(viewer, &format));
730: if (format == PETSC_VIEWER_ASCII_INFO) {
731: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO run parameters:\n"));
732: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO phase: %d \n", mat_mkl_cpardiso->phase));
733: for (PetscInt i = 1; i <= 64; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO iparm[%d]: %d \n", i, mat_mkl_cpardiso->iparm[i - 1]));
734: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO maxfct: %d \n", mat_mkl_cpardiso->maxfct));
735: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO mnum: %d \n", mat_mkl_cpardiso->mnum));
736: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO mtype: %d \n", mat_mkl_cpardiso->mtype));
737: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO n: %d \n", mat_mkl_cpardiso->n));
738: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO nrhs: %d \n", mat_mkl_cpardiso->nrhs));
739: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO msglvl: %d \n", mat_mkl_cpardiso->msglvl));
740: }
741: }
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: static PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info)
746: {
747: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
749: PetscFunctionBegin;
750: info->block_size = 1.0;
751: info->nz_allocated = mat_mkl_cpardiso->nz + 0.0;
752: info->nz_unneeded = 0.0;
753: info->assemblies = 0.0;
754: info->mallocs = 0.0;
755: info->memory = 0.0;
756: info->fill_ratio_given = 0;
757: info->fill_ratio_needed = 0;
758: info->factor_mallocs = 0;
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: static PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F, PetscInt icntl, PetscInt ival)
763: {
764: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
766: PetscFunctionBegin;
767: if (icntl <= 64) {
768: mat_mkl_cpardiso->iparm[icntl - 1] = ival;
769: } else {
770: if (icntl == 65) mkl_set_num_threads((int)ival);
771: else if (icntl == 66) mat_mkl_cpardiso->maxfct = ival;
772: else if (icntl == 67) mat_mkl_cpardiso->mnum = ival;
773: else if (icntl == 68) mat_mkl_cpardiso->msglvl = ival;
774: else if (icntl == 69) mat_mkl_cpardiso->mtype = ival;
775: }
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
779: /*@
780: MatMkl_CPardisoSetCntl - Set MKL Cluster PARDISO parameters
781: <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html>
783: Logically Collective
785: Input Parameters:
786: + F - the factored matrix obtained by calling `MatGetFactor()`
787: . icntl - index of MKL Cluster PARDISO parameter
788: - ival - value of MKL Cluster PARDISO parameter
790: Options Database Key:
791: . -mat_mkl_cpardiso_icntl ival - set the option numbered icntl to ival
793: Level: intermediate
795: Note:
796: This routine cannot be used if you are solving the linear system with `TS`, `SNES`, or `KSP`, only if you directly call `MatGetFactor()` so use the options
797: database approach when working with `TS`, `SNES`, or `KSP`. See `MATSOLVERMKL_CPARDISO` for the options
799: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MATMPIAIJ`, `MATSOLVERMKL_CPARDISO`
800: @*/
801: PetscErrorCode MatMkl_CPardisoSetCntl(Mat F, PetscInt icntl, PetscInt ival)
802: {
803: PetscFunctionBegin;
804: PetscTryMethod(F, "MatMkl_CPardisoSetCntl_C", (Mat, PetscInt, PetscInt), (F, icntl, ival));
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: /*MC
809: MATSOLVERMKL_CPARDISO - A matrix type providing direct solvers (LU) for parallel matrices via the external package MKL Cluster PARDISO
810: <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html>
812: Works with `MATMPIAIJ` matrices
814: Use `-pc_type lu` `-pc_factor_mat_solver_type mkl_cpardiso` to use this direct solver
816: Options Database Keys:
817: + -mat_mkl_cpardiso_65 - Suggested number of threads to use within MKL Cluster PARDISO
818: . -mat_mkl_cpardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
819: . -mat_mkl_cpardiso_67 - Indicates the actual matrix for the solution phase
820: . -mat_mkl_cpardiso_68 - Message level information, use 1 to get detailed information on the solver options
821: . -mat_mkl_cpardiso_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
822: . -mat_mkl_cpardiso_1 - Use default values
823: . -mat_mkl_cpardiso_2 - Fill-in reducing ordering for the input matrix
824: . -mat_mkl_cpardiso_4 - Preconditioned CGS/CG
825: . -mat_mkl_cpardiso_5 - User permutation
826: . -mat_mkl_cpardiso_6 - Write solution on x
827: . -mat_mkl_cpardiso_8 - Iterative refinement step
828: . -mat_mkl_cpardiso_10 - Pivoting perturbation
829: . -mat_mkl_cpardiso_11 - Scaling vectors
830: . -mat_mkl_cpardiso_12 - Solve with transposed or conjugate transposed matrix A
831: . -mat_mkl_cpardiso_13 - Improved accuracy using (non-) symmetric weighted matching
832: . -mat_mkl_cpardiso_18 - Numbers of non-zero elements
833: . -mat_mkl_cpardiso_19 - Report number of floating point operations
834: . -mat_mkl_cpardiso_21 - Pivoting for symmetric indefinite matrices
835: . -mat_mkl_cpardiso_24 - Parallel factorization control
836: . -mat_mkl_cpardiso_25 - Parallel forward/backward solve control
837: . -mat_mkl_cpardiso_27 - Matrix checker
838: . -mat_mkl_cpardiso_31 - Partial solve and computing selected components of the solution vectors
839: . -mat_mkl_cpardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
840: - -mat_mkl_cpardiso_60 - Intel MKL Cluster PARDISO mode
842: Level: beginner
844: Notes:
845: Use `-mat_mkl_cpardiso_68 1` to display the number of threads the solver is using. MKL does not provide a way to directly access this
846: information.
848: For more information on the options check
849: <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html>
851: .seealso: [](ch_matrices), `Mat`, `PCFactorSetMatSolverType()`, `MatSolverType`, `MatMkl_CPardisoSetCntl()`, `MatGetFactor()`, `MATSOLVERMKL_PARDISO`
852: M*/
854: static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type)
855: {
856: PetscFunctionBegin;
857: *type = MATSOLVERMKL_CPARDISO;
858: PetscFunctionReturn(PETSC_SUCCESS);
859: }
861: /* MatGetFactor for MPI AIJ matrices */
862: static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A, MatFactorType ftype, Mat *F)
863: {
864: Mat B;
865: Mat_MKL_CPARDISO *mat_mkl_cpardiso;
866: PetscBool isSeqAIJ, isMPIBAIJ, isMPISBAIJ;
868: PetscFunctionBegin;
869: /* Create the factorization matrix */
871: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSeqAIJ));
872: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIBAIJ, &isMPIBAIJ));
873: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISBAIJ, &isMPISBAIJ));
874: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
875: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
876: PetscCall(PetscStrallocpy("mkl_cpardiso", &((PetscObject)B)->type_name));
877: PetscCall(MatSetUp(B));
879: PetscCall(PetscNew(&mat_mkl_cpardiso));
881: if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO;
882: else if (isMPIBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO;
883: else if (isMPISBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO;
884: else mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO;
886: if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO;
887: else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_CPARDISO;
888: B->ops->destroy = MatDestroy_MKL_CPARDISO;
890: B->ops->view = MatView_MKL_CPARDISO;
891: B->ops->getinfo = MatGetInfo_MKL_CPARDISO;
893: B->factortype = ftype;
894: B->assembled = PETSC_TRUE; /* required by -ksp_view */
896: B->data = mat_mkl_cpardiso;
898: /* set solvertype */
899: PetscCall(PetscFree(B->solvertype));
900: PetscCall(PetscStrallocpy(MATSOLVERMKL_CPARDISO, &B->solvertype));
902: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_mkl_cpardiso));
903: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMkl_CPardisoSetCntl_C", MatMkl_CPardisoSetCntl_MKL_CPARDISO));
904: PetscCall(PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso));
906: *F = B;
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void)
911: {
912: PetscFunctionBegin;
913: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
914: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
915: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPIBAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
916: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPISBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_mpiaij_mkl_cpardiso));
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }