1 /*
2 Provides an interface to the PaStiX sparse solver
3 */
4 #include <../src/mat/impls/aij/seq/aij.h>
5 #include <../src/mat/impls/aij/mpi/mpiaij.h>
6 #include <../src/mat/impls/sbaij/seq/sbaij.h>
7 #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
8
9 #if defined(PETSC_USE_COMPLEX)
10 #define _H_COMPLEX
11 #endif
12
13 #include <pastix.h>
14
15 #if defined(PETSC_USE_COMPLEX)
16 #if defined(PETSC_USE_REAL_SINGLE)
17 #define SPM_FLTTYPE SpmComplex32
18 #else
19 #define SPM_FLTTYPE SpmComplex64
20 #endif
21 #else /* PETSC_USE_COMPLEX */
22
23 #if defined(PETSC_USE_REAL_SINGLE)
24 #define SPM_FLTTYPE SpmFloat
25 #else
26 #define SPM_FLTTYPE SpmDouble
27 #endif
28
29 #endif /* PETSC_USE_COMPLEX */
30
31 typedef PetscScalar PastixScalar;
32
33 typedef struct Mat_Pastix_ {
34 pastix_data_t *pastix_data; /* Pastix data storage structure */
35 MPI_Comm comm; /* MPI Communicator used to initialize pastix */
36 spmatrix_t *spm; /* SPM matrix structure */
37 MatStructure matstruc; /* DIFFERENT_NONZERO_PATTERN if uninitialized, SAME otherwise */
38 PetscScalar *rhs; /* Right-hand-side member */
39 PetscInt rhsnbr; /* Right-hand-side number */
40 pastix_int_t iparm[IPARM_SIZE]; /* Integer parameters */
41 double dparm[DPARM_SIZE]; /* Floating point parameters */
42 } Mat_Pastix;
43
44 /*
45 convert PETSc matrix to SPM structure
46
47 input:
48 A - matrix in aij, baij or sbaij format
49 reuse - MAT_INITIAL_MATRIX: spaces are allocated and values are set for the triple
50 MAT_REUSE_MATRIX: only the values in v array are updated
51 output:
52 spm - The SPM built from A
53 */
MatConvertToSPM(Mat A,MatReuse reuse,Mat_Pastix * pastix)54 static PetscErrorCode MatConvertToSPM(Mat A, MatReuse reuse, Mat_Pastix *pastix)
55 {
56 Mat A_loc, A_aij;
57 const PetscInt *row, *col;
58 PetscInt n, i;
59 const PetscScalar *val;
60 PetscBool ismpiaij, isseqaij, ismpisbaij, isseqsbaij;
61 PetscBool flag;
62 spmatrix_t spm2, *spm = NULL;
63 int spm_err;
64
65 PetscFunctionBegin;
66 /* Get A datas */
67 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
68 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
69 /* TODO: Block Aij should be handled with dof in spm */
70 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQSBAIJ, &isseqsbaij));
71 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPISBAIJ, &ismpisbaij));
72
73 if (isseqsbaij || ismpisbaij) PetscCall(MatConvert(A, MATAIJ, reuse, &A_aij));
74 else A_aij = A;
75
76 if (ismpiaij || ismpisbaij) PetscCall(MatMPIAIJGetLocalMat(A_aij, MAT_INITIAL_MATRIX, &A_loc));
77 else if (isseqaij || isseqsbaij) A_loc = A_aij;
78 else SETERRQ(PetscObjectComm((PetscObject)A_aij), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
79
80 /* Use getRowIJ and the trick CSC/CSR instead of GetColumnIJ for performance */
81 PetscCall(MatGetRowIJ(A_loc, 0, PETSC_FALSE, PETSC_FALSE, &n, &row, &col, &flag));
82 PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_SUP, "GetRowIJ failed");
83 PetscCall(MatSeqAIJGetArrayRead(A_loc, &val));
84
85 PetscCall(PetscMalloc1(1, &spm));
86 PetscStackCallExternalVoid("spmInitDist", spmInitDist(spm, pastix->comm));
87
88 spm->n = n;
89 spm->nnz = row[n];
90 spm->fmttype = SpmCSR;
91 spm->flttype = SPM_FLTTYPE;
92 spm->replicated = !(A->rmap->n != A->rmap->N);
93
94 PetscStackCallExternalVoid("spmUpdateComputedFields", spmUpdateComputedFields(spm));
95 PetscStackCallExternalVoid("spmAlloc", spmAlloc(spm));
96
97 /* Get data distribution */
98 if (!spm->replicated) {
99 for (i = A->rmap->rstart; i < A->rmap->rend; i++) spm->loc2glob[i - A->rmap->rstart] = i;
100 }
101
102 /* Copy arrays */
103 PetscCall(PetscArraycpy(spm->colptr, col, spm->nnz));
104 PetscCall(PetscArraycpy(spm->rowptr, row, spm->n + 1));
105 PetscCall(PetscArraycpy((PetscScalar *)spm->values, val, spm->nnzexp));
106 PetscCall(MatRestoreRowIJ(A_loc, 0, PETSC_FALSE, PETSC_FALSE, &n, &row, &col, &flag));
107 PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_SUP, "RestoreRowIJ failed");
108 PetscCall(MatSeqAIJRestoreArrayRead(A_loc, &val));
109 if (ismpiaij || ismpisbaij) PetscCall(MatDestroy(&A_loc));
110
111 /*
112 PaStiX works only with CSC matrices for now, so let's lure him by telling him
113 that the PETSc CSR matrix is a CSC matrix.
114 Note that this is not available yet for Hermitian matrices and LL^h/LDL^h factorizations.
115 */
116 {
117 spm_int_t *tmp;
118 spm->fmttype = SpmCSC;
119 tmp = spm->colptr;
120 spm->colptr = spm->rowptr;
121 spm->rowptr = tmp;
122 pastix->iparm[IPARM_TRANSPOSE_SOLVE] = PastixTrans;
123 }
124
125 /* Update matrix to be in PaStiX format */
126 PetscStackCallExternalVoid("spmCheckAndCorrect", spm_err = spmCheckAndCorrect(spm, &spm2));
127 if (spm_err != 0) {
128 PetscStackCallExternalVoid("spmExit", spmExit(spm));
129 *spm = spm2;
130 }
131
132 if (isseqsbaij || ismpisbaij) PetscCall(MatDestroy(&A_aij));
133
134 pastix->spm = spm;
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
138 /*
139 Call clean step of PaStiX if initialized
140 Free the SpM matrix and the PaStiX instance.
141 */
MatDestroy_PaStiX(Mat A)142 static PetscErrorCode MatDestroy_PaStiX(Mat A)
143 {
144 Mat_Pastix *pastix = (Mat_Pastix *)A->data;
145
146 PetscFunctionBegin;
147 /* Finalize SPM (matrix handler of PaStiX) */
148 if (pastix->spm) {
149 PetscStackCallExternalVoid("spmExit", spmExit(pastix->spm));
150 PetscCall(PetscFree(pastix->spm));
151 }
152
153 /* clear composed functions */
154 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
155
156 /* Finalize PaStiX */
157 if (pastix->pastix_data) pastixFinalize(&pastix->pastix_data);
158
159 /* Deallocate PaStiX structure */
160 PetscCall(PetscFree(A->data));
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
164 /*
165 Gather right-hand side.
166 Call for Solve step.
167 Scatter solution.
168 */
MatSolve_PaStiX(Mat A,Vec b,Vec x)169 static PetscErrorCode MatSolve_PaStiX(Mat A, Vec b, Vec x)
170 {
171 Mat_Pastix *pastix = (Mat_Pastix *)A->data;
172 const PetscScalar *bptr;
173 PetscInt ldrhs;
174
175 PetscFunctionBegin;
176 pastix->rhsnbr = 1;
177 ldrhs = pastix->spm->n;
178
179 PetscCall(VecCopy(b, x));
180 PetscCall(VecGetArray(x, &pastix->rhs));
181 PetscCall(VecGetArrayRead(b, &bptr));
182
183 /* solve phase */
184 /*-------------*/
185 PetscCheck(pastix->pastix_data, PETSC_COMM_SELF, PETSC_ERR_SUP, "PaStiX hasn't been initialized");
186 PetscCallExternal(pastix_task_solve, pastix->pastix_data, ldrhs, pastix->rhsnbr, pastix->rhs, ldrhs);
187 PetscCallExternal(pastix_task_refine, pastix->pastix_data, ldrhs, pastix->rhsnbr, (PetscScalar *)bptr, ldrhs, pastix->rhs, ldrhs);
188
189 PetscCall(VecRestoreArray(x, &pastix->rhs));
190 PetscCall(VecRestoreArrayRead(b, &bptr));
191 PetscFunctionReturn(PETSC_SUCCESS);
192 }
193
194 /*
195 Numeric factorisation using PaStiX solver.
196
197 input:
198 F - PETSc matrix that contains PaStiX interface.
199 A - PETSc matrix in aij, bail or sbaij format
200 */
MatFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo * info)201 static PetscErrorCode MatFactorNumeric_PaStiX(Mat F, Mat A, const MatFactorInfo *info)
202 {
203 Mat_Pastix *pastix = (Mat_Pastix *)F->data;
204
205 PetscFunctionBegin;
206 F->ops->solve = MatSolve_PaStiX;
207
208 /* Perform Numerical Factorization */
209 PetscCheck(pastix->pastix_data, PETSC_COMM_SELF, PETSC_ERR_SUP, "PaStiX hasn't been initialized");
210 PetscCallExternal(pastix_task_numfact, pastix->pastix_data, pastix->spm);
211
212 F->assembled = PETSC_TRUE;
213 PetscFunctionReturn(PETSC_SUCCESS);
214 }
215
MatLUFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo * info)216 static PetscErrorCode MatLUFactorNumeric_PaStiX(Mat F, Mat A, const MatFactorInfo *info)
217 {
218 Mat_Pastix *pastix = (Mat_Pastix *)F->data;
219
220 PetscFunctionBegin;
221 PetscCheck(pastix->iparm[IPARM_FACTORIZATION] == PastixFactGETRF, PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Incorrect factorization type for symbolic and numerical factorization by PaStiX");
222 pastix->iparm[IPARM_FACTORIZATION] = PastixFactGETRF;
223 PetscCall(MatFactorNumeric_PaStiX(F, A, info));
224 PetscFunctionReturn(PETSC_SUCCESS);
225 }
226
MatCholeskyFactorNumeric_PaStiX(Mat F,Mat A,const MatFactorInfo * info)227 static PetscErrorCode MatCholeskyFactorNumeric_PaStiX(Mat F, Mat A, const MatFactorInfo *info)
228 {
229 Mat_Pastix *pastix = (Mat_Pastix *)F->data;
230
231 PetscFunctionBegin;
232 PetscCheck(pastix->iparm[IPARM_FACTORIZATION] == PastixFactSYTRF, PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Incorrect factorization type for symbolic and numerical factorization by PaStiX");
233 pastix->iparm[IPARM_FACTORIZATION] = PastixFactSYTRF;
234 PetscCall(MatFactorNumeric_PaStiX(F, A, info));
235 PetscFunctionReturn(PETSC_SUCCESS);
236 }
237
238 /*
239 Perform Ordering step and Symbolic Factorization step
240
241 Note the PETSc r and c permutations are ignored
242 input:
243 F - PETSc matrix that contains PaStiX interface.
244 A - matrix in aij, bail or sbaij format
245 r - permutation ?
246 c - TODO
247 info - Information about the factorization to perform.
248 output:
249 pastix_data - This instance will be updated with the SolverMatrix allocated.
250 */
MatFactorSymbolic_PaStiX(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)251 static PetscErrorCode MatFactorSymbolic_PaStiX(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
252 {
253 Mat_Pastix *pastix = (Mat_Pastix *)F->data;
254
255 PetscFunctionBegin;
256 pastix->matstruc = DIFFERENT_NONZERO_PATTERN;
257
258 /* Initialise SPM structure */
259 PetscCall(MatConvertToSPM(A, MAT_INITIAL_MATRIX, pastix));
260
261 /* Ordering - Symbolic factorization - Build SolverMatrix */
262 PetscCheck(pastix->pastix_data, PETSC_COMM_SELF, PETSC_ERR_SUP, "PaStiX hasn't been initialized");
263 PetscCallExternal(pastix_task_analyze, pastix->pastix_data, pastix->spm);
264 PetscFunctionReturn(PETSC_SUCCESS);
265 }
266
MatLUFactorSymbolic_PaStiX(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)267 static PetscErrorCode MatLUFactorSymbolic_PaStiX(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
268 {
269 Mat_Pastix *pastix = (Mat_Pastix *)F->data;
270
271 PetscFunctionBegin;
272 pastix->iparm[IPARM_FACTORIZATION] = PastixFactGETRF;
273 PetscCall(MatFactorSymbolic_PaStiX(F, A, r, c, info));
274 PetscFunctionReturn(PETSC_SUCCESS);
275 }
276
277 /* Note the PETSc r permutation is ignored */
MatCholeskyFactorSymbolic_PaStiX(Mat F,Mat A,IS r,const MatFactorInfo * info)278 static PetscErrorCode MatCholeskyFactorSymbolic_PaStiX(Mat F, Mat A, IS r, const MatFactorInfo *info)
279 {
280 Mat_Pastix *pastix = (Mat_Pastix *)F->data;
281
282 PetscFunctionBegin;
283 /* Warning: Cholesky in PETSc wrapper does not handle (complex) Hermitian matrices.
284 The factorization type can be forced using the parameter
285 mat_pastix_factorization (see enum pastix_factotype_e in
286 https://solverstack.gitlabpages.inria.fr/pastix/group__pastix__api.html). */
287 pastix->iparm[IPARM_FACTORIZATION] = PastixFactSYTRF;
288 PetscCall(MatFactorSymbolic_PaStiX(F, A, r, NULL, info));
289 PetscFunctionReturn(PETSC_SUCCESS);
290 }
291
MatView_PaStiX(Mat A,PetscViewer viewer)292 static PetscErrorCode MatView_PaStiX(Mat A, PetscViewer viewer)
293 {
294 PetscBool isascii;
295 PetscViewerFormat format;
296
297 PetscFunctionBegin;
298 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
299 if (isascii) {
300 PetscCall(PetscViewerGetFormat(viewer, &format));
301 if (format == PETSC_VIEWER_ASCII_INFO) {
302 Mat_Pastix *pastix = (Mat_Pastix *)A->data;
303 spmatrix_t *spm = pastix->spm;
304 PetscCheck(!spm, PETSC_COMM_SELF, PETSC_ERR_SUP, "Sparse matrix isn't initialized");
305
306 PetscCall(PetscViewerASCIIPrintf(viewer, "PaStiX run parameters:\n"));
307 PetscCall(PetscViewerASCIIPrintf(viewer, " Matrix type : %s \n", ((spm->mtxtype == SpmSymmetric) ? "Symmetric" : "Unsymmetric")));
308 PetscCall(PetscViewerASCIIPrintf(viewer, " Level of printing (0,1,2): %ld \n", (long)pastix->iparm[IPARM_VERBOSE]));
309 PetscCall(PetscViewerASCIIPrintf(viewer, " Number of refinements iterations : %ld \n", (long)pastix->iparm[IPARM_NBITER]));
310 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error : %e \n", pastix->dparm[DPARM_RELATIVE_ERROR]));
311 if (pastix->iparm[IPARM_VERBOSE] > 0) spmPrintInfo(spm, stdout);
312 }
313 }
314 PetscFunctionReturn(PETSC_SUCCESS);
315 }
316
317 /*MC
318 MATSOLVERPASTIX - A solver package providing direct solvers (LU) for distributed
319 and sequential matrices via the external package PaStiX.
320
321 Use `./configure --download-hwloc --download-metis --download-ptscotch --download-pastix --download-netlib-lapack [or MKL for ex. --with-blaslapack-dir=${MKLROOT}]`
322 to have PETSc installed with PasTiX.
323
324 Use `-pc_type lu` `-pc_factor_mat_solver_type pastix` to use this direct solver.
325
326 Options Database Keys:
327 -mat_pastix_verbose <0,1,2> - print level of information messages from PaStiX
328 -mat_pastix_factorization <0,1,2,3,4> - Factorization algorithm (Cholesky, LDL^t, LU, LL^t, LDL^h)
329 -mat_pastix_itermax <integer> - Maximum number of iterations during refinement
330 -mat_pastix_epsilon_refinement <double> - Epsilon for refinement
331 -mat_pastix_epsilon_magn_ctrl <double> - Epsilon for magnitude control
332 -mat_pastix_ordering <0,1> - Ordering (Scotch or Metis)
333 -mat_pastix_thread_nbr <integer> - Set the numbers of threads for each MPI process
334 -mat_pastix_scheduler <0,1,2,3,4> - Scheduler (sequential, static, parsec, starpu, dynamic)
335 -mat_pastix_compress_when <0,1,2,3> - When to compress (never, minimal-theory, just-in-time, supernodes)
336 -mat_pastix_compress_tolerance <double> - Tolerance for low-rank kernels.
337
338 Notes:
339 This only works for matrices with symmetric nonzero structure, if you pass it a matrix with
340 nonsymmetric structure PasTiX, and hence, PETSc return with an error.
341
342 Level: beginner
343
344 .seealso: [](ch_matrices), `Mat`, `PCFactorSetMatSolverType()`, `MatSolverType`, `MatGetFactor()`
345 M*/
346
MatGetInfo_PaStiX(Mat A,MatInfoType flag,MatInfo * info)347 static PetscErrorCode MatGetInfo_PaStiX(Mat A, MatInfoType flag, MatInfo *info)
348 {
349 Mat_Pastix *pastix = (Mat_Pastix *)A->data;
350
351 PetscFunctionBegin;
352 info->block_size = 1.0;
353 info->nz_allocated = pastix->iparm[IPARM_ALLOCATED_TERMS];
354 info->nz_used = pastix->iparm[IPARM_NNZEROS];
355 info->nz_unneeded = 0.0;
356 info->assemblies = 0.0;
357 info->mallocs = 0.0;
358 info->memory = 0.0;
359 info->fill_ratio_given = 0;
360 info->fill_ratio_needed = 0;
361 info->factor_mallocs = 0;
362 PetscFunctionReturn(PETSC_SUCCESS);
363 }
364
MatFactorGetSolverType_PaStiX(Mat A,MatSolverType * type)365 static PetscErrorCode MatFactorGetSolverType_PaStiX(Mat A, MatSolverType *type)
366 {
367 PetscFunctionBegin;
368 *type = MATSOLVERPASTIX;
369 PetscFunctionReturn(PETSC_SUCCESS);
370 }
371
372 /* Sets PaStiX options from the options database */
MatSetFromOptions_PaStiX(Mat A)373 static PetscErrorCode MatSetFromOptions_PaStiX(Mat A)
374 {
375 Mat_Pastix *pastix = (Mat_Pastix *)A->data;
376 pastix_int_t *iparm = pastix->iparm;
377 double *dparm = pastix->dparm;
378 PetscInt icntl;
379 PetscReal dcntl;
380 PetscBool set;
381
382 PetscFunctionBegin;
383 PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "PaStiX Options", "Mat");
384
385 iparm[IPARM_VERBOSE] = 0;
386 iparm[IPARM_ITERMAX] = 20;
387
388 PetscCall(PetscOptionsRangeInt("-mat_pastix_verbose", "iparm[IPARM_VERBOSE] : level of printing (0 to 2)", "None", iparm[IPARM_VERBOSE], &icntl, &set, 0, 2));
389 if (set) iparm[IPARM_VERBOSE] = (pastix_int_t)icntl;
390
391 PetscCall(PetscOptionsRangeInt("-mat_pastix_factorization", "iparm[IPARM_FACTORIZATION]: Factorization algorithm", "None", iparm[IPARM_FACTORIZATION], &icntl, &set, 0, 4));
392 if (set) iparm[IPARM_FACTORIZATION] = (pastix_int_t)icntl;
393
394 PetscCall(PetscOptionsBoundedInt("-mat_pastix_itermax", "iparm[IPARM_ITERMAX]: Max iterations", "None", iparm[IPARM_ITERMAX], &icntl, &set, 1));
395 if (set) iparm[IPARM_ITERMAX] = (pastix_int_t)icntl;
396
397 PetscCall(PetscOptionsBoundedReal("-mat_pastix_epsilon_refinement", "dparm[DPARM_EPSILON_REFINEMENT]: Epsilon refinement", "None", dparm[DPARM_EPSILON_REFINEMENT], &dcntl, &set, -1.));
398 if (set) dparm[DPARM_EPSILON_REFINEMENT] = (double)dcntl;
399
400 PetscCall(PetscOptionsReal("-mat_pastix_epsilon_magn_ctrl", "dparm[DPARM_EPSILON_MAGN_CTRL]: Epsilon magnitude control", "None", dparm[DPARM_EPSILON_MAGN_CTRL], &dcntl, &set));
401 if (set) dparm[DPARM_EPSILON_MAGN_CTRL] = (double)dcntl;
402
403 PetscCall(PetscOptionsRangeInt("-mat_pastix_ordering", "iparm[IPARM_ORDERING]: Ordering algorithm", "None", iparm[IPARM_ORDERING], &icntl, &set, 0, 2));
404 if (set) iparm[IPARM_ORDERING] = (pastix_int_t)icntl;
405
406 PetscCall(PetscOptionsBoundedInt("-mat_pastix_thread_nbr", "iparm[IPARM_THREAD_NBR]: Number of thread by MPI node", "None", iparm[IPARM_THREAD_NBR], &icntl, &set, -1));
407 if (set) iparm[IPARM_THREAD_NBR] = (pastix_int_t)icntl;
408
409 PetscCall(PetscOptionsRangeInt("-mat_pastix_scheduler", "iparm[IPARM_SCHEDULER]: Scheduler", "None", iparm[IPARM_SCHEDULER], &icntl, &set, 0, 4));
410 if (set) iparm[IPARM_SCHEDULER] = (pastix_int_t)icntl;
411
412 PetscCall(PetscOptionsRangeInt("-mat_pastix_compress_when", "iparm[IPARM_COMPRESS_WHEN]: When to compress", "None", iparm[IPARM_COMPRESS_WHEN], &icntl, &set, 0, 3));
413 if (set) iparm[IPARM_COMPRESS_WHEN] = (pastix_int_t)icntl;
414
415 PetscCall(PetscOptionsBoundedReal("-mat_pastix_compress_tolerance", "dparm[DPARM_COMPRESS_TOLERANCE]: Tolerance for low-rank kernels", "None", dparm[DPARM_COMPRESS_TOLERANCE], &dcntl, &set, 0.));
416 if (set) dparm[DPARM_COMPRESS_TOLERANCE] = (double)dcntl;
417
418 PetscOptionsEnd();
419 PetscFunctionReturn(PETSC_SUCCESS);
420 }
421
MatGetFactor_pastix(Mat A,MatFactorType ftype,Mat * F,const char * mattype)422 static PetscErrorCode MatGetFactor_pastix(Mat A, MatFactorType ftype, Mat *F, const char *mattype)
423 {
424 Mat B;
425 Mat_Pastix *pastix;
426
427 PetscFunctionBegin;
428 /* Create the factorization matrix */
429 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
430 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
431 PetscCall(PetscStrallocpy(MATSOLVERPASTIX, &((PetscObject)B)->type_name));
432 PetscCall(MatSetUp(B));
433
434 PetscCheck(ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_CHOLESKY, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported by PaStiX");
435
436 /* set solvertype */
437 PetscCall(PetscFree(B->solvertype));
438 PetscCall(PetscStrallocpy(MATSOLVERPASTIX, &B->solvertype));
439
440 B->ops->lufactorsymbolic = MatLUFactorSymbolic_PaStiX;
441 B->ops->lufactornumeric = MatLUFactorNumeric_PaStiX;
442 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_PaStiX;
443 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_PaStiX;
444 B->ops->view = MatView_PaStiX;
445 B->ops->getinfo = MatGetInfo_PaStiX;
446 B->ops->destroy = MatDestroy_PaStiX;
447
448 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_PaStiX));
449
450 B->factortype = ftype;
451
452 /* Create the pastix structure */
453 PetscCall(PetscNew(&pastix));
454 B->data = (void *)pastix;
455
456 /* Call to set default pastix options */
457 PetscStackCallExternalVoid("pastixInitParam", pastixInitParam(pastix->iparm, pastix->dparm));
458 PetscCall(MatSetFromOptions_PaStiX(B));
459
460 /* Get PETSc communicator */
461 PetscCall(PetscObjectGetComm((PetscObject)A, &pastix->comm));
462
463 /* Initialise PaStiX structure */
464 pastix->iparm[IPARM_SCOTCH_MT] = 0;
465 PetscStackCallExternalVoid("pastixInit", pastixInit(&pastix->pastix_data, pastix->comm, pastix->iparm, pastix->dparm));
466
467 /* Warning: Cholesky in PETSc wrapper does not handle (complex) Hermitian matrices.
468 The factorization type can be forced using the parameter
469 mat_pastix_factorization (see enum pastix_factotype_e in
470 https://solverstack.gitlabpages.inria.fr/pastix/group__pastix__api.html). */
471 pastix->iparm[IPARM_FACTORIZATION] = ftype == MAT_FACTOR_CHOLESKY ? PastixFactSYTRF : PastixFactGETRF;
472
473 *F = B;
474 PetscFunctionReturn(PETSC_SUCCESS);
475 }
476
MatGetFactor_mpiaij_pastix(Mat A,MatFactorType ftype,Mat * F)477 static PetscErrorCode MatGetFactor_mpiaij_pastix(Mat A, MatFactorType ftype, Mat *F)
478 {
479 PetscFunctionBegin;
480 PetscCheck(ftype == MAT_FACTOR_LU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
481 PetscCall(MatGetFactor_pastix(A, ftype, F, MATMPIAIJ));
482 PetscFunctionReturn(PETSC_SUCCESS);
483 }
484
MatGetFactor_seqaij_pastix(Mat A,MatFactorType ftype,Mat * F)485 static PetscErrorCode MatGetFactor_seqaij_pastix(Mat A, MatFactorType ftype, Mat *F)
486 {
487 PetscFunctionBegin;
488 PetscCheck(ftype == MAT_FACTOR_LU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc AIJ matrices with PaStiX Cholesky, use SBAIJ matrix");
489 PetscCall(MatGetFactor_pastix(A, ftype, F, MATSEQAIJ));
490 PetscFunctionReturn(PETSC_SUCCESS);
491 }
492
MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat * F)493 static PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A, MatFactorType ftype, Mat *F)
494 {
495 PetscFunctionBegin;
496 PetscCheck(ftype == MAT_FACTOR_CHOLESKY, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
497 PetscCall(MatGetFactor_pastix(A, ftype, F, MATMPISBAIJ));
498 PetscFunctionReturn(PETSC_SUCCESS);
499 }
500
MatGetFactor_seqsbaij_pastix(Mat A,MatFactorType ftype,Mat * F)501 static PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat A, MatFactorType ftype, Mat *F)
502 {
503 PetscFunctionBegin;
504 PetscCheck(ftype == MAT_FACTOR_CHOLESKY, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");
505 PetscCall(MatGetFactor_pastix(A, ftype, F, MATSEQSBAIJ));
506 PetscFunctionReturn(PETSC_SUCCESS);
507 }
508
MatSolverTypeRegister_Pastix(void)509 PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Pastix(void)
510 {
511 PetscFunctionBegin;
512 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_pastix));
513 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_pastix));
514 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATMPISBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_mpisbaij_pastix));
515 PetscCall(MatSolverTypeRegister(MATSOLVERPASTIX, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_pastix));
516 PetscFunctionReturn(PETSC_SUCCESS);
517 }
518