Actual source code: mkl_cpardiso.c
2: #include <petscsys.h>
3: #include <../src/mat/impls/aij/mpi/mpiaij.h>
4: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
6: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
7: #define MKL_ILP64
8: #endif
9: #include <mkl.h>
10: #include <mkl_cluster_sparse_solver.h>
12: /*
13: * Possible mkl_cpardiso phases that controls the execution of the solver.
14: * For more information check mkl_cpardiso manual.
15: */
16: #define JOB_ANALYSIS 11
17: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
18: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
19: #define JOB_NUMERICAL_FACTORIZATION 22
20: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
21: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
22: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
23: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
24: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
25: #define JOB_RELEASE_OF_LU_MEMORY 0
26: #define JOB_RELEASE_OF_ALL_MEMORY -1
28: #define IPARM_SIZE 64
29: #define INT_TYPE MKL_INT
31: static const char *Err_MSG_CPardiso(int errNo)
32: {
33: switch (errNo) {
34: case -1:
35: return "input inconsistent";
36: break;
37: case -2:
38: return "not enough memory";
39: break;
40: case -3:
41: return "reordering problem";
42: break;
43: case -4:
44: return "zero pivot, numerical factorization or iterative refinement problem";
45: break;
46: case -5:
47: return "unclassified (internal) error";
48: break;
49: case -6:
50: return "preordering failed (matrix types 11, 13 only)";
51: break;
52: case -7:
53: return "diagonal matrix problem";
54: break;
55: case -8:
56: return "32-bit integer overflow problem";
57: break;
58: case -9:
59: return "not enough memory for OOC";
60: break;
61: case -10:
62: return "problems with opening OOC temporary files";
63: break;
64: case -11:
65: return "read/write problems with the OOC data file";
66: break;
67: default:
68: return "unknown error";
69: }
70: }
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 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 sparce 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: 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: 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;
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: PetscErrorCode MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
215: {
216: const PetscInt *ai, *aj, *bi, *bj, *garray, bs = A->rmap->bs, bs2 = bs * bs, m = A->rmap->n / bs, *ajj, *bjj;
217: PetscInt rstart, nz, i, j, countA, countB;
218: PetscInt *row, *col;
219: const PetscScalar *av, *bv;
220: PetscScalar *val;
221: Mat_MPIBAIJ *mat = (Mat_MPIBAIJ *)A->data;
222: Mat_SeqBAIJ *aa = (Mat_SeqBAIJ *)(mat->A)->data;
223: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ *)(mat->B)->data;
224: PetscInt colA_start, jB, jcol;
226: PetscFunctionBegin;
227: ai = aa->i;
228: aj = aa->j;
229: bi = bb->i;
230: bj = bb->j;
231: rstart = A->rmap->rstart / bs;
232: av = aa->a;
233: bv = bb->a;
235: garray = mat->garray;
237: if (reuse == MAT_INITIAL_MATRIX) {
238: nz = aa->nz + bb->nz;
239: *nnz = nz;
240: PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz * bs2, &val));
241: *r = row;
242: *c = col;
243: *v = val;
244: } else {
245: row = *r;
246: col = *c;
247: val = *v;
248: }
250: nz = 0;
251: for (i = 0; i < m; i++) {
252: row[i] = nz + 1;
253: countA = ai[i + 1] - ai[i];
254: countB = bi[i + 1] - bi[i];
255: ajj = aj + ai[i]; /* ptr to the beginning of this row */
256: bjj = bj + bi[i];
258: /* B part, smaller col index */
259: colA_start = rstart + (countA > 0 ? ajj[0] : 0); /* the smallest global col index of A */
260: jB = 0;
261: for (j = 0; j < countB; j++) {
262: jcol = garray[bjj[j]];
263: if (jcol > colA_start) break;
264: col[nz++] = jcol + 1;
265: }
266: jB = j;
267: PetscCall(PetscArraycpy(val, bv, jB * bs2));
268: val += jB * bs2;
269: bv += jB * bs2;
271: /* A part */
272: for (j = 0; j < countA; j++) col[nz++] = rstart + ajj[j] + 1;
273: PetscCall(PetscArraycpy(val, av, countA * bs2));
274: val += countA * bs2;
275: av += countA * bs2;
277: /* B part, larger col index */
278: for (j = jB; j < countB; j++) col[nz++] = garray[bjj[j]] + 1;
279: PetscCall(PetscArraycpy(val, bv, (countB - jB) * bs2));
280: val += (countB - jB) * bs2;
281: bv += (countB - jB) * bs2;
282: }
283: row[m] = nz + 1;
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
289: {
290: const PetscInt *ai, *aj, *bi, *bj, *garray, bs = A->rmap->bs, bs2 = bs * bs, m = A->rmap->n / bs, *ajj, *bjj;
291: PetscInt rstart, nz, i, j, countA, countB;
292: PetscInt *row, *col;
293: const PetscScalar *av, *bv;
294: PetscScalar *val;
295: Mat_MPISBAIJ *mat = (Mat_MPISBAIJ *)A->data;
296: Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ *)(mat->A)->data;
297: Mat_SeqBAIJ *bb = (Mat_SeqBAIJ *)(mat->B)->data;
299: PetscFunctionBegin;
300: ai = aa->i;
301: aj = aa->j;
302: bi = bb->i;
303: bj = bb->j;
304: rstart = A->rmap->rstart / bs;
305: av = aa->a;
306: bv = bb->a;
308: garray = mat->garray;
310: if (reuse == MAT_INITIAL_MATRIX) {
311: nz = aa->nz + bb->nz;
312: *nnz = nz;
313: PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz * bs2, &val));
314: *r = row;
315: *c = col;
316: *v = val;
317: } else {
318: row = *r;
319: col = *c;
320: val = *v;
321: }
323: nz = 0;
324: for (i = 0; i < m; i++) {
325: row[i] = nz + 1;
326: countA = ai[i + 1] - ai[i];
327: countB = bi[i + 1] - bi[i];
328: ajj = aj + ai[i]; /* ptr to the beginning of this row */
329: bjj = bj + bi[i];
331: /* A part */
332: for (j = 0; j < countA; j++) col[nz++] = rstart + ajj[j] + 1;
333: PetscCall(PetscArraycpy(val, av, countA * bs2));
334: val += countA * bs2;
335: av += countA * bs2;
337: /* B part, larger col index */
338: for (j = 0; j < countB; j++) col[nz++] = garray[bjj[j]] + 1;
339: PetscCall(PetscArraycpy(val, bv, countB * bs2));
340: val += countB * bs2;
341: bv += countB * bs2;
342: }
343: row[m] = nz + 1;
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*
349: * Free memory for Mat_MKL_CPARDISO structure and pointers to objects.
350: */
351: PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A)
352: {
353: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
354: MPI_Comm comm;
356: PetscFunctionBegin;
357: /* Terminate instance, deallocate memories */
358: if (mat_mkl_cpardiso->CleanUp) {
359: mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;
361: 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,
362: mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err);
363: }
365: if (mat_mkl_cpardiso->ConvertToTriples != MatCopy_seqaij_seqaij_MKL_CPARDISO) PetscCall(PetscFree3(mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja, mat_mkl_cpardiso->a));
366: comm = MPI_Comm_f2c(mat_mkl_cpardiso->comm_mkl_cpardiso);
367: PetscCallMPI(MPI_Comm_free(&comm));
368: PetscCall(PetscFree(A->data));
370: /* clear composed functions */
371: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
372: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMkl_CPardisoSetCntl_C", NULL));
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: /*
377: * Computes Ax = b
378: */
379: PetscErrorCode MatSolve_MKL_CPARDISO(Mat A, Vec b, Vec x)
380: {
381: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)(A)->data;
382: PetscScalar *xarray;
383: const PetscScalar *barray;
385: PetscFunctionBegin;
386: mat_mkl_cpardiso->nrhs = 1;
387: PetscCall(VecGetArray(x, &xarray));
388: PetscCall(VecGetArrayRead(b, &barray));
390: /* solve phase */
391: mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
392: 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,
393: 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);
395: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
397: PetscCall(VecRestoreArray(x, &xarray));
398: PetscCall(VecRestoreArrayRead(b, &barray));
399: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A, Vec b, Vec x)
404: {
405: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
407: PetscFunctionBegin;
408: #if defined(PETSC_USE_COMPLEX)
409: mat_mkl_cpardiso->iparm[12 - 1] = 1;
410: #else
411: mat_mkl_cpardiso->iparm[12 - 1] = 2;
412: #endif
413: PetscCall(MatSolve_MKL_CPARDISO(A, b, x));
414: mat_mkl_cpardiso->iparm[12 - 1] = 0;
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A, Mat B, Mat X)
419: {
420: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)(A)->data;
421: PetscScalar *xarray;
422: const PetscScalar *barray;
424: PetscFunctionBegin;
425: PetscCall(MatGetSize(B, NULL, (PetscInt *)&mat_mkl_cpardiso->nrhs));
427: if (mat_mkl_cpardiso->nrhs > 0) {
428: PetscCall(MatDenseGetArrayRead(B, &barray));
429: PetscCall(MatDenseGetArray(X, &xarray));
431: PetscCheck(barray != xarray, PETSC_COMM_SELF, PETSC_ERR_SUP, "B and X cannot share the same memory location");
433: /* solve phase */
434: mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
435: 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,
436: 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);
437: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
438: PetscCall(MatDenseRestoreArrayRead(B, &barray));
439: PetscCall(MatDenseRestoreArray(X, &xarray));
440: }
441: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: /*
446: * LU Decomposition
447: */
448: PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F, Mat A, const MatFactorInfo *info)
449: {
450: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)(F)->data;
452: PetscFunctionBegin;
453: mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
454: 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));
456: mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION;
457: 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,
458: 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);
459: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
461: mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
462: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: /* Sets mkl_cpardiso options from the options database */
467: PetscErrorCode MatSetFromOptions_MKL_CPARDISO(Mat F, Mat A)
468: {
469: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
470: PetscInt icntl, threads;
471: PetscBool flg;
473: PetscFunctionBegin;
474: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "MKL_CPARDISO Options", "Mat");
475: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_65", "Suggested number of threads to use within MKL_CPARDISO", "None", threads, &threads, &flg));
476: if (flg) mkl_set_num_threads((int)threads);
478: 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));
479: if (flg) mat_mkl_cpardiso->maxfct = icntl;
481: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_67", "Indicates the actual matrix for the solution phase", "None", mat_mkl_cpardiso->mnum, &icntl, &flg));
482: if (flg) mat_mkl_cpardiso->mnum = icntl;
484: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_68", "Message level information", "None", mat_mkl_cpardiso->msglvl, &icntl, &flg));
485: if (flg) mat_mkl_cpardiso->msglvl = icntl;
487: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_69", "Defines the matrix type", "None", mat_mkl_cpardiso->mtype, &icntl, &flg));
488: if (flg) mat_mkl_cpardiso->mtype = icntl;
489: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_1", "Use default values", "None", mat_mkl_cpardiso->iparm[0], &icntl, &flg));
491: if (flg && icntl != 0) {
492: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_2", "Fill-in reducing ordering for the input matrix", "None", mat_mkl_cpardiso->iparm[1], &icntl, &flg));
493: if (flg) mat_mkl_cpardiso->iparm[1] = icntl;
495: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_4", "Preconditioned CGS/CG", "None", mat_mkl_cpardiso->iparm[3], &icntl, &flg));
496: if (flg) mat_mkl_cpardiso->iparm[3] = icntl;
498: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_5", "User permutation", "None", mat_mkl_cpardiso->iparm[4], &icntl, &flg));
499: if (flg) mat_mkl_cpardiso->iparm[4] = icntl;
501: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_6", "Write solution on x", "None", mat_mkl_cpardiso->iparm[5], &icntl, &flg));
502: if (flg) mat_mkl_cpardiso->iparm[5] = icntl;
504: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_8", "Iterative refinement step", "None", mat_mkl_cpardiso->iparm[7], &icntl, &flg));
505: if (flg) mat_mkl_cpardiso->iparm[7] = icntl;
507: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_10", "Pivoting perturbation", "None", mat_mkl_cpardiso->iparm[9], &icntl, &flg));
508: if (flg) mat_mkl_cpardiso->iparm[9] = icntl;
510: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_11", "Scaling vectors", "None", mat_mkl_cpardiso->iparm[10], &icntl, &flg));
511: if (flg) mat_mkl_cpardiso->iparm[10] = icntl;
513: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_12", "Solve with transposed or conjugate transposed matrix A", "None", mat_mkl_cpardiso->iparm[11], &icntl, &flg));
514: if (flg) mat_mkl_cpardiso->iparm[11] = icntl;
516: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_13", "Improved accuracy using (non-) symmetric weighted matching", "None", mat_mkl_cpardiso->iparm[12], &icntl, &flg));
517: if (flg) mat_mkl_cpardiso->iparm[12] = icntl;
519: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_18", "Numbers of non-zero elements", "None", mat_mkl_cpardiso->iparm[17], &icntl, &flg));
520: if (flg) mat_mkl_cpardiso->iparm[17] = icntl;
522: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_19", "Report number of floating point operations", "None", mat_mkl_cpardiso->iparm[18], &icntl, &flg));
523: if (flg) mat_mkl_cpardiso->iparm[18] = icntl;
525: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_21", "Pivoting for symmetric indefinite matrices", "None", mat_mkl_cpardiso->iparm[20], &icntl, &flg));
526: if (flg) mat_mkl_cpardiso->iparm[20] = icntl;
528: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_24", "Parallel factorization control", "None", mat_mkl_cpardiso->iparm[23], &icntl, &flg));
529: if (flg) mat_mkl_cpardiso->iparm[23] = icntl;
531: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_25", "Parallel forward/backward solve control", "None", mat_mkl_cpardiso->iparm[24], &icntl, &flg));
532: if (flg) mat_mkl_cpardiso->iparm[24] = icntl;
534: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_27", "Matrix checker", "None", mat_mkl_cpardiso->iparm[26], &icntl, &flg));
535: if (flg) mat_mkl_cpardiso->iparm[26] = icntl;
537: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_31", "Partial solve and computing selected components of the solution vectors", "None", mat_mkl_cpardiso->iparm[30], &icntl, &flg));
538: if (flg) mat_mkl_cpardiso->iparm[30] = icntl;
540: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_34", "Optimal number of threads for conditional numerical reproducibility (CNR) mode", "None", mat_mkl_cpardiso->iparm[33], &icntl, &flg));
541: if (flg) mat_mkl_cpardiso->iparm[33] = icntl;
543: PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_60", "Intel MKL_CPARDISO mode", "None", mat_mkl_cpardiso->iparm[59], &icntl, &flg));
544: if (flg) mat_mkl_cpardiso->iparm[59] = icntl;
545: }
547: PetscOptionsEnd();
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso)
552: {
553: PetscInt bs;
554: PetscBool match;
555: PetscMPIInt size;
556: MPI_Comm comm;
558: PetscFunctionBegin;
560: PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)A), &comm));
561: PetscCallMPI(MPI_Comm_size(comm, &size));
562: mat_mkl_cpardiso->comm_mkl_cpardiso = MPI_Comm_c2f(comm);
564: mat_mkl_cpardiso->CleanUp = PETSC_FALSE;
565: mat_mkl_cpardiso->maxfct = 1;
566: mat_mkl_cpardiso->mnum = 1;
567: mat_mkl_cpardiso->n = A->rmap->N;
568: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
569: mat_mkl_cpardiso->msglvl = 0;
570: mat_mkl_cpardiso->nrhs = 1;
571: mat_mkl_cpardiso->err = 0;
572: mat_mkl_cpardiso->phase = -1;
573: #if defined(PETSC_USE_COMPLEX)
574: mat_mkl_cpardiso->mtype = 13;
575: #else
576: mat_mkl_cpardiso->mtype = 11;
577: #endif
579: #if defined(PETSC_USE_REAL_SINGLE)
580: mat_mkl_cpardiso->iparm[27] = 1;
581: #else
582: mat_mkl_cpardiso->iparm[27] = 0;
583: #endif
585: mat_mkl_cpardiso->iparm[0] = 1; /* Solver default parameters overridden with provided by iparm */
586: mat_mkl_cpardiso->iparm[1] = 2; /* Use METIS for fill-in reordering */
587: mat_mkl_cpardiso->iparm[5] = 0; /* Write solution into x */
588: mat_mkl_cpardiso->iparm[7] = 2; /* Max number of iterative refinement steps */
589: mat_mkl_cpardiso->iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
590: mat_mkl_cpardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
591: mat_mkl_cpardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
592: mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
593: mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
594: mat_mkl_cpardiso->iparm[26] = 1; /* Check input data for correctness */
596: mat_mkl_cpardiso->iparm[39] = 0;
597: if (size > 1) {
598: mat_mkl_cpardiso->iparm[39] = 2;
599: mat_mkl_cpardiso->iparm[40] = A->rmap->rstart;
600: mat_mkl_cpardiso->iparm[41] = A->rmap->rend - 1;
601: }
602: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &match, MATMPIBAIJ, MATMPISBAIJ, ""));
603: if (match) {
604: PetscCall(MatGetBlockSize(A, &bs));
605: mat_mkl_cpardiso->iparm[36] = bs;
606: mat_mkl_cpardiso->iparm[40] /= bs;
607: mat_mkl_cpardiso->iparm[41] /= bs;
608: mat_mkl_cpardiso->iparm[40]++;
609: mat_mkl_cpardiso->iparm[41]++;
610: mat_mkl_cpardiso->iparm[34] = 0; /* Fortran style */
611: } else {
612: mat_mkl_cpardiso->iparm[34] = 1; /* C style */
613: }
615: mat_mkl_cpardiso->perm = 0;
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*
620: * Symbolic decomposition. Mkl_Pardiso analysis phase.
621: */
622: PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
623: {
624: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
626: PetscFunctionBegin;
627: mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
629: /* Set MKL_CPARDISO options from the options database */
630: PetscCall(MatSetFromOptions_MKL_CPARDISO(F, A));
631: 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));
633: mat_mkl_cpardiso->n = A->rmap->N;
634: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
636: /* analysis phase */
637: mat_mkl_cpardiso->phase = JOB_ANALYSIS;
639: 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,
640: 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);
642: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL_CPARDISO: err=%d, msg = \"%s\".Check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
644: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
645: F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO;
646: F->ops->solve = MatSolve_MKL_CPARDISO;
647: F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO;
648: F->ops->matsolve = MatMatSolve_MKL_CPARDISO;
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_CPARDISO(Mat F, Mat A, IS perm, const MatFactorInfo *info)
653: {
654: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
656: PetscFunctionBegin;
657: mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;
659: /* Set MKL_CPARDISO options from the options database */
660: PetscCall(MatSetFromOptions_MKL_CPARDISO(F, A));
661: 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));
663: mat_mkl_cpardiso->n = A->rmap->N;
664: if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
665: PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead");
666: if (A->spd == PETSC_BOOL3_TRUE) mat_mkl_cpardiso->mtype = 2;
667: else mat_mkl_cpardiso->mtype = -2;
669: /* analysis phase */
670: mat_mkl_cpardiso->phase = JOB_ANALYSIS;
672: 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,
673: 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);
675: PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL_CPARDISO: err=%d, msg = \"%s\".Check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
677: mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
678: F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_CPARDISO;
679: F->ops->solve = MatSolve_MKL_CPARDISO;
680: F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO;
681: F->ops->matsolve = MatMatSolve_MKL_CPARDISO;
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer)
686: {
687: PetscBool iascii;
688: PetscViewerFormat format;
689: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
690: PetscInt i;
692: PetscFunctionBegin;
693: /* check if matrix is mkl_cpardiso type */
694: if (A->ops->solve != MatSolve_MKL_CPARDISO) PetscFunctionReturn(PETSC_SUCCESS);
696: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
697: if (iascii) {
698: PetscCall(PetscViewerGetFormat(viewer, &format));
699: if (format == PETSC_VIEWER_ASCII_INFO) {
700: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO run parameters:\n"));
701: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO phase: %d \n", mat_mkl_cpardiso->phase));
702: for (i = 1; i <= 64; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO iparm[%d]: %d \n", i, mat_mkl_cpardiso->iparm[i - 1]));
703: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO maxfct: %d \n", mat_mkl_cpardiso->maxfct));
704: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO mnum: %d \n", mat_mkl_cpardiso->mnum));
705: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO mtype: %d \n", mat_mkl_cpardiso->mtype));
706: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO n: %d \n", mat_mkl_cpardiso->n));
707: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO nrhs: %d \n", mat_mkl_cpardiso->nrhs));
708: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL_CPARDISO msglvl: %d \n", mat_mkl_cpardiso->msglvl));
709: }
710: }
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info)
715: {
716: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
718: PetscFunctionBegin;
719: info->block_size = 1.0;
720: info->nz_allocated = mat_mkl_cpardiso->nz + 0.0;
721: info->nz_unneeded = 0.0;
722: info->assemblies = 0.0;
723: info->mallocs = 0.0;
724: info->memory = 0.0;
725: info->fill_ratio_given = 0;
726: info->fill_ratio_needed = 0;
727: info->factor_mallocs = 0;
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F, PetscInt icntl, PetscInt ival)
732: {
733: Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
735: PetscFunctionBegin;
736: if (icntl <= 64) {
737: mat_mkl_cpardiso->iparm[icntl - 1] = ival;
738: } else {
739: if (icntl == 65) mkl_set_num_threads((int)ival);
740: else if (icntl == 66) mat_mkl_cpardiso->maxfct = ival;
741: else if (icntl == 67) mat_mkl_cpardiso->mnum = ival;
742: else if (icntl == 68) mat_mkl_cpardiso->msglvl = ival;
743: else if (icntl == 69) mat_mkl_cpardiso->mtype = ival;
744: }
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: /*@
749: MatMkl_CPardisoSetCntl - Set Mkl_Pardiso parameters
751: Logically Collective
753: Input Parameters:
754: + F - the factored matrix obtained by calling `MatGetFactor()`
755: . icntl - index of Mkl_Pardiso parameter
756: - ival - value of Mkl_Pardiso parameter
758: Options Database Key:
759: . -mat_mkl_cpardiso_<icntl> <ival> - set the option numbered icntl to ival
761: Level: Intermediate
763: Note:
764: 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
765: database approach when working with `TS`, `SNES`, or `KSP`. See `MATSOLVERMKL_CPARDISO` for the options
767: References:
768: . * - Mkl_Pardiso Users' Guide
770: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MATMPIAIJ`, `MATSOLVERMKL_CPARDISO`
771: @*/
772: PetscErrorCode MatMkl_CPardisoSetCntl(Mat F, PetscInt icntl, PetscInt ival)
773: {
774: PetscFunctionBegin;
775: PetscTryMethod(F, "MatMkl_CPardisoSetCntl_C", (Mat, PetscInt, PetscInt), (F, icntl, ival));
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
779: /*MC
780: MATSOLVERMKL_CPARDISO - A matrix type providing direct solvers (LU) for parallel matrices via the external package MKL_CPARDISO.
782: Works with `MATMPIAIJ` matrices
784: Use `-pc_type lu` `-pc_factor_mat_solver_type mkl_cpardiso` to use this direct solver
786: Options Database Keys:
787: + -mat_mkl_cpardiso_65 - Suggested number of threads to use within MKL_CPARDISO
788: . -mat_mkl_cpardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
789: . -mat_mkl_cpardiso_67 - Indicates the actual matrix for the solution phase
790: . -mat_mkl_cpardiso_68 - Message level information, use 1 to get detailed information on the solver options
791: . -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
792: . -mat_mkl_cpardiso_1 - Use default values
793: . -mat_mkl_cpardiso_2 - Fill-in reducing ordering for the input matrix
794: . -mat_mkl_cpardiso_4 - Preconditioned CGS/CG
795: . -mat_mkl_cpardiso_5 - User permutation
796: . -mat_mkl_cpardiso_6 - Write solution on x
797: . -mat_mkl_cpardiso_8 - Iterative refinement step
798: . -mat_mkl_cpardiso_10 - Pivoting perturbation
799: . -mat_mkl_cpardiso_11 - Scaling vectors
800: . -mat_mkl_cpardiso_12 - Solve with transposed or conjugate transposed matrix A
801: . -mat_mkl_cpardiso_13 - Improved accuracy using (non-) symmetric weighted matching
802: . -mat_mkl_cpardiso_18 - Numbers of non-zero elements
803: . -mat_mkl_cpardiso_19 - Report number of floating point operations
804: . -mat_mkl_cpardiso_21 - Pivoting for symmetric indefinite matrices
805: . -mat_mkl_cpardiso_24 - Parallel factorization control
806: . -mat_mkl_cpardiso_25 - Parallel forward/backward solve control
807: . -mat_mkl_cpardiso_27 - Matrix checker
808: . -mat_mkl_cpardiso_31 - Partial solve and computing selected components of the solution vectors
809: . -mat_mkl_cpardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
810: - -mat_mkl_cpardiso_60 - Intel MKL_CPARDISO mode
812: Level: beginner
814: Notes:
815: 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
816: information.
818: For more information on the options check the MKL_CPARDISO manual
820: .seealso: [](ch_matrices), `Mat`, `PCFactorSetMatSolverType()`, `MatSolverType`, `MatMkl_CPardisoSetCntl()`, `MatGetFactor()`
821: M*/
823: static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type)
824: {
825: PetscFunctionBegin;
826: *type = MATSOLVERMKL_CPARDISO;
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: /* MatGetFactor for MPI AIJ matrices */
831: static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A, MatFactorType ftype, Mat *F)
832: {
833: Mat B;
834: Mat_MKL_CPARDISO *mat_mkl_cpardiso;
835: PetscBool isSeqAIJ, isMPIBAIJ, isMPISBAIJ;
837: PetscFunctionBegin;
838: /* Create the factorization matrix */
840: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSeqAIJ));
841: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIBAIJ, &isMPIBAIJ));
842: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISBAIJ, &isMPISBAIJ));
843: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
844: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
845: PetscCall(PetscStrallocpy("mkl_cpardiso", &((PetscObject)B)->type_name));
846: PetscCall(MatSetUp(B));
848: PetscCall(PetscNew(&mat_mkl_cpardiso));
850: if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO;
851: else if (isMPIBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO;
852: else if (isMPISBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO;
853: else mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO;
855: if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO;
856: else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_CPARDISO;
857: B->ops->destroy = MatDestroy_MKL_CPARDISO;
859: B->ops->view = MatView_MKL_CPARDISO;
860: B->ops->getinfo = MatGetInfo_MKL_CPARDISO;
862: B->factortype = ftype;
863: B->assembled = PETSC_TRUE; /* required by -ksp_view */
865: B->data = mat_mkl_cpardiso;
867: /* set solvertype */
868: PetscCall(PetscFree(B->solvertype));
869: PetscCall(PetscStrallocpy(MATSOLVERMKL_CPARDISO, &B->solvertype));
871: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_mkl_cpardiso));
872: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMkl_CPardisoSetCntl_C", MatMkl_CPardisoSetCntl_MKL_CPARDISO));
873: PetscCall(PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso));
875: *F = B;
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void)
880: {
881: PetscFunctionBegin;
882: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
883: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
884: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPIBAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
885: PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPISBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_mpiaij_mkl_cpardiso));
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }