1 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
2 #include <../src/mat/impls/aij/mpi/mpiaij.h>
3 #include <StrumpackSparseSolver.h>
4
MatGetDiagonal_STRUMPACK(Mat A,Vec v)5 static PetscErrorCode MatGetDiagonal_STRUMPACK(Mat A,Vec v)
6 {
7 PetscFunctionBegin;
8 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK factor");
9 PetscFunctionReturn(0);
10 }
11
MatDestroy_STRUMPACK(Mat A)12 static PetscErrorCode MatDestroy_STRUMPACK(Mat A)
13 {
14 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
15 PetscErrorCode ierr;
16 PetscBool flg;
17
18 PetscFunctionBegin;
19 /* Deallocate STRUMPACK storage */
20 PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(S));
21 ierr = PetscFree(A->spptr);CHKERRQ(ierr);
22 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr);
23 if (flg) {
24 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr);
25 } else {
26 ierr = MatDestroy_MPIAIJ(A);CHKERRQ(ierr);
27 }
28
29 /* clear composed functions */
30 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);CHKERRQ(ierr);
31 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetReordering_C",NULL);CHKERRQ(ierr);
32 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetColPerm_C",NULL);CHKERRQ(ierr);
33 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSRelTol_C",NULL);CHKERRQ(ierr);
34 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSAbsTol_C",NULL);CHKERRQ(ierr);
35 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMaxRank_C",NULL);CHKERRQ(ierr);
36 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSLeafSize_C",NULL);CHKERRQ(ierr);
37 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMinSepSize_C",NULL);CHKERRQ(ierr);
38
39 PetscFunctionReturn(0);
40 }
41
42
MatSTRUMPACKSetReordering_STRUMPACK(Mat F,MatSTRUMPACKReordering reordering)43 static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F,MatSTRUMPACKReordering reordering)
44 {
45 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
46
47 PetscFunctionBegin;
48 PetscStackCall("STRUMPACK_reordering_method",STRUMPACK_set_reordering_method(*S,(STRUMPACK_REORDERING_STRATEGY)reordering));
49 PetscFunctionReturn(0);
50 }
51
52 /*@
53 MatSTRUMPACKSetReordering - Set STRUMPACK fill-reducing reordering
54
55 Input Parameters:
56 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
57 - reordering - the code to be used to find the fill-reducing reordering
58 Possible values: NATURAL=0 METIS=1 PARMETIS=2 SCOTCH=3 PTSCOTCH=4 RCM=5
59
60 Options Database:
61 . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
62
63 Level: beginner
64
65 References:
66 . STRUMPACK manual
67
68 .seealso: MatGetFactor()
69 @*/
MatSTRUMPACKSetReordering(Mat F,MatSTRUMPACKReordering reordering)70 PetscErrorCode MatSTRUMPACKSetReordering(Mat F,MatSTRUMPACKReordering reordering)
71 {
72 PetscErrorCode ierr;
73
74 PetscFunctionBegin;
75 PetscValidHeaderSpecific(F,MAT_CLASSID,1);
76 PetscValidLogicalCollectiveEnum(F,reordering,2);
77 ierr = PetscTryMethod(F,"MatSTRUMPACKSetReordering_C",(Mat,MatSTRUMPACKReordering),(F,reordering));CHKERRQ(ierr);
78 PetscFunctionReturn(0);
79 }
80
MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)81 static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)
82 {
83 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
84
85 PetscFunctionBegin;
86 PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,cperm ? 5 : 0));
87 PetscFunctionReturn(0);
88 }
89
90 /*@
91 MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
92 Logically Collective on Mat
93
94 Input Parameters:
95 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
96 - cperm - whether or not to permute (internally) the columns of the matrix
97
98 Options Database:
99 . -mat_strumpack_colperm <cperm>
100
101 Level: beginner
102
103 References:
104 . STRUMPACK manual
105
106 .seealso: MatGetFactor()
107 @*/
MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)108 PetscErrorCode MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)
109 {
110 PetscErrorCode ierr;
111
112 PetscFunctionBegin;
113 PetscValidHeaderSpecific(F,MAT_CLASSID,1);
114 PetscValidLogicalCollectiveBool(F,cperm,2);
115 ierr = PetscTryMethod(F,"MatSTRUMPACKSetColPerm_C",(Mat,PetscBool),(F,cperm));CHKERRQ(ierr);
116 PetscFunctionReturn(0);
117 }
118
MatSTRUMPACKSetHSSRelTol_STRUMPACK(Mat F,PetscReal rtol)119 static PetscErrorCode MatSTRUMPACKSetHSSRelTol_STRUMPACK(Mat F,PetscReal rtol)
120 {
121 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
122
123 PetscFunctionBegin;
124 PetscStackCall("STRUMPACK_set_HSS_rel_tol", STRUMPACK_set_HSS_rel_tol(*S,rtol));
125 PetscFunctionReturn(0);
126 }
127
128 /*@
129 MatSTRUMPACKSetHSSRelTol - Set STRUMPACK relative tolerance for HSS compression
130 Logically Collective on Mat
131
132 Input Parameters:
133 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
134 - rtol - relative compression tolerance
135
136 Options Database:
137 . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
138
139 Level: beginner
140
141 References:
142 . STRUMPACK manual
143
144 .seealso: MatGetFactor()
145 @*/
MatSTRUMPACKSetHSSRelTol(Mat F,PetscReal rtol)146 PetscErrorCode MatSTRUMPACKSetHSSRelTol(Mat F,PetscReal rtol)
147 {
148 PetscErrorCode ierr;
149
150 PetscFunctionBegin;
151 PetscValidHeaderSpecific(F,MAT_CLASSID,1);
152 PetscValidLogicalCollectiveReal(F,rtol,2);
153 ierr = PetscTryMethod(F,"MatSTRUMPACKSetHSSRelTol_C",(Mat,PetscReal),(F,rtol));CHKERRQ(ierr);
154 PetscFunctionReturn(0);
155 }
156
MatSTRUMPACKSetHSSAbsTol_STRUMPACK(Mat F,PetscReal atol)157 static PetscErrorCode MatSTRUMPACKSetHSSAbsTol_STRUMPACK(Mat F,PetscReal atol)
158 {
159 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
160
161 PetscFunctionBegin;
162 PetscStackCall("STRUMPACK_set_HSS_abs_tol", STRUMPACK_set_HSS_abs_tol(*S,atol));
163 PetscFunctionReturn(0);
164 }
165
166 /*@
167 MatSTRUMPACKSetHSSAbsTol - Set STRUMPACK absolute tolerance for HSS compression
168 Logically Collective on Mat
169
170 Input Parameters:
171 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
172 - atol - absolute compression tolerance
173
174 Options Database:
175 . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
176
177 Level: beginner
178
179 References:
180 . STRUMPACK manual
181
182 .seealso: MatGetFactor()
183 @*/
MatSTRUMPACKSetHSSAbsTol(Mat F,PetscReal atol)184 PetscErrorCode MatSTRUMPACKSetHSSAbsTol(Mat F,PetscReal atol)
185 {
186 PetscErrorCode ierr;
187
188 PetscFunctionBegin;
189 PetscValidHeaderSpecific(F,MAT_CLASSID,1);
190 PetscValidLogicalCollectiveReal(F,atol,2);
191 ierr = PetscTryMethod(F,"MatSTRUMPACKSetHSSAbsTol_C",(Mat,PetscReal),(F,atol));CHKERRQ(ierr);
192 PetscFunctionReturn(0);
193 }
194
MatSTRUMPACKSetHSSMaxRank_STRUMPACK(Mat F,PetscInt hssmaxrank)195 static PetscErrorCode MatSTRUMPACKSetHSSMaxRank_STRUMPACK(Mat F,PetscInt hssmaxrank)
196 {
197 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
198
199 PetscFunctionBegin;
200 PetscStackCall("STRUMPACK_set_HSS_max_rank", STRUMPACK_set_HSS_max_rank(*S,hssmaxrank));
201 PetscFunctionReturn(0);
202 }
203
204 /*@
205 MatSTRUMPACKSetHSSMaxRank - Set STRUMPACK maximum HSS rank
206 Logically Collective on Mat
207
208 Input Parameters:
209 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
210 - hssmaxrank - maximum rank used in low-rank approximation
211
212 Options Database:
213 . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
214
215 Level: beginner
216
217 References:
218 . STRUMPACK manual
219
220 .seealso: MatGetFactor()
221 @*/
MatSTRUMPACKSetHSSMaxRank(Mat F,PetscInt hssmaxrank)222 PetscErrorCode MatSTRUMPACKSetHSSMaxRank(Mat F,PetscInt hssmaxrank)
223 {
224 PetscErrorCode ierr;
225
226 PetscFunctionBegin;
227 PetscValidHeaderSpecific(F,MAT_CLASSID,1);
228 PetscValidLogicalCollectiveInt(F,hssmaxrank,2);
229 ierr = PetscTryMethod(F,"MatSTRUMPACKSetHSSMaxRank_C",(Mat,PetscInt),(F,hssmaxrank));CHKERRQ(ierr);
230 PetscFunctionReturn(0);
231 }
232
MatSTRUMPACKSetHSSLeafSize_STRUMPACK(Mat F,PetscInt leaf_size)233 static PetscErrorCode MatSTRUMPACKSetHSSLeafSize_STRUMPACK(Mat F,PetscInt leaf_size)
234 {
235 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
236
237 PetscFunctionBegin;
238 PetscStackCall("STRUMPACK_set_HSS_leaf_size", STRUMPACK_set_HSS_leaf_size(*S,leaf_size));
239 PetscFunctionReturn(0);
240 }
241
242 /*@
243 MatSTRUMPACKSetHSSLeafSize - Set STRUMPACK HSS leaf size
244 Logically Collective on Mat
245
246 Input Parameters:
247 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
248 - leaf_size - Size of diagonal blocks in HSS approximation
249
250 Options Database:
251 . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
252
253 Level: beginner
254
255 References:
256 . STRUMPACK manual
257
258 .seealso: MatGetFactor()
259 @*/
MatSTRUMPACKSetHSSLeafSize(Mat F,PetscInt leaf_size)260 PetscErrorCode MatSTRUMPACKSetHSSLeafSize(Mat F,PetscInt leaf_size)
261 {
262 PetscErrorCode ierr;
263
264 PetscFunctionBegin;
265 PetscValidHeaderSpecific(F,MAT_CLASSID,1);
266 PetscValidLogicalCollectiveInt(F,leaf_size,2);
267 ierr = PetscTryMethod(F,"MatSTRUMPACKSetHSSLeafSize_C",(Mat,PetscInt),(F,leaf_size));CHKERRQ(ierr);
268 PetscFunctionReturn(0);
269 }
270
MatSTRUMPACKSetHSSMinSepSize_STRUMPACK(Mat F,PetscInt hssminsize)271 static PetscErrorCode MatSTRUMPACKSetHSSMinSepSize_STRUMPACK(Mat F,PetscInt hssminsize)
272 {
273 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
274
275 PetscFunctionBegin;
276 PetscStackCall("STRUMPACK_set_HSS_min_sep_size", STRUMPACK_set_HSS_min_sep_size(*S,hssminsize));
277 PetscFunctionReturn(0);
278 }
279
280 /*@
281 MatSTRUMPACKSetHSSMinSepSize - Set STRUMPACK minimum separator size for low-rank approximation
282 Logically Collective on Mat
283
284 Input Parameters:
285 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
286 - hssminsize - minimum dense matrix size for low-rank approximation
287
288 Options Database:
289 . -mat_strumpack_hss_min_sep_size <hssminsize>
290
291 Level: beginner
292
293 References:
294 . STRUMPACK manual
295
296 .seealso: MatGetFactor()
297 @*/
MatSTRUMPACKSetHSSMinSepSize(Mat F,PetscInt hssminsize)298 PetscErrorCode MatSTRUMPACKSetHSSMinSepSize(Mat F,PetscInt hssminsize)
299 {
300 PetscErrorCode ierr;
301
302 PetscFunctionBegin;
303 PetscValidHeaderSpecific(F,MAT_CLASSID,1);
304 PetscValidLogicalCollectiveInt(F,hssminsize,2);
305 ierr = PetscTryMethod(F,"MatSTRUMPACKSetHSSMinSepSize_C",(Mat,PetscInt),(F,hssminsize));CHKERRQ(ierr);
306 PetscFunctionReturn(0);
307 }
308
MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)309 static PetscErrorCode MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)
310 {
311 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
312 STRUMPACK_RETURN_CODE sp_err;
313 PetscErrorCode ierr;
314 const PetscScalar *bptr;
315 PetscScalar *xptr;
316
317 PetscFunctionBegin;
318 ierr = VecGetArray(x,&xptr);CHKERRQ(ierr);
319 ierr = VecGetArrayRead(b_mpi,&bptr);CHKERRQ(ierr);
320
321 PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(*S,(PetscScalar*)bptr,xptr,0));
322 switch (sp_err) {
323 case STRUMPACK_SUCCESS: break;
324 case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
325 case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
326 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
327 }
328 ierr = VecRestoreArray(x,&xptr);CHKERRQ(ierr);
329 ierr = VecRestoreArrayRead(b_mpi,&bptr);CHKERRQ(ierr);
330 PetscFunctionReturn(0);
331 }
332
MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)333 static PetscErrorCode MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)
334 {
335 PetscErrorCode ierr;
336 PetscBool flg;
337
338 PetscFunctionBegin;
339 ierr = PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr);
340 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
341 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr);
342 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
343 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK() is not implemented yet");
344 PetscFunctionReturn(0);
345 }
346
MatView_Info_STRUMPACK(Mat A,PetscViewer viewer)347 static PetscErrorCode MatView_Info_STRUMPACK(Mat A,PetscViewer viewer)
348 {
349 PetscErrorCode ierr;
350
351 PetscFunctionBegin;
352 /* check if matrix is strumpack type */
353 if (A->ops->solve != MatSolve_STRUMPACK) PetscFunctionReturn(0);
354 ierr = PetscViewerASCIIPrintf(viewer,"STRUMPACK sparse solver!\n");CHKERRQ(ierr);
355 PetscFunctionReturn(0);
356 }
357
MatView_STRUMPACK(Mat A,PetscViewer viewer)358 static PetscErrorCode MatView_STRUMPACK(Mat A,PetscViewer viewer)
359 {
360 PetscErrorCode ierr;
361 PetscBool iascii;
362 PetscViewerFormat format;
363
364 PetscFunctionBegin;
365 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
366 if (iascii) {
367 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
368 if (format == PETSC_VIEWER_ASCII_INFO) {
369 ierr = MatView_Info_STRUMPACK(A,viewer);CHKERRQ(ierr);
370 }
371 }
372 PetscFunctionReturn(0);
373 }
374
MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo * info)375 static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo *info)
376 {
377 STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
378 STRUMPACK_RETURN_CODE sp_err;
379 Mat_SeqAIJ *A_d,*A_o;
380 Mat_MPIAIJ *mat;
381 PetscErrorCode ierr;
382 PetscInt M=A->rmap->N,m=A->rmap->n;
383 PetscBool flg;
384
385 PetscFunctionBegin;
386 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
387 if (flg) { /* A is MATMPIAIJ */
388 mat = (Mat_MPIAIJ*)A->data;
389 A_d = (Mat_SeqAIJ*)(mat->A)->data;
390 A_o = (Mat_SeqAIJ*)(mat->B)->data;
391 PetscStackCall("STRUMPACK_set_MPIAIJ_matrix",STRUMPACK_set_MPIAIJ_matrix(*S,&m,A_d->i,A_d->j,A_d->a,A_o->i,A_o->j,A_o->a,mat->garray));
392 } else { /* A is MATSEQAIJ */
393 A_d = (Mat_SeqAIJ*)A->data;
394 PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(*S,&M,A_d->i,A_d->j,A_d->a,0));
395 }
396
397 /* Reorder and Factor the matrix. */
398 /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
399 PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(*S));
400 PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(*S));
401 switch (sp_err) {
402 case STRUMPACK_SUCCESS: break;
403 case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
404 case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
405 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
406 }
407 F->assembled = PETSC_TRUE;
408 F->preallocated = PETSC_TRUE;
409 PetscFunctionReturn(0);
410 }
411
MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)412 static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
413 {
414 PetscFunctionBegin;
415 F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
416 F->ops->solve = MatSolve_STRUMPACK;
417 F->ops->matsolve = MatMatSolve_STRUMPACK;
418 PetscFunctionReturn(0);
419 }
420
MatFactorGetSolverType_aij_strumpack(Mat A,MatSolverType * type)421 static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A,MatSolverType *type)
422 {
423 PetscFunctionBegin;
424 *type = MATSOLVERSTRUMPACK;
425 PetscFunctionReturn(0);
426 }
427
428 /*MC
429 MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
430 and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.
431
432 Consult the STRUMPACK-sparse manual for more info.
433
434 Use
435 ./configure --download-strumpack
436 to have PETSc installed with STRUMPACK
437
438 Use
439 -pc_type lu -pc_factor_mat_solver_type strumpack
440 to use this as an exact (direct) solver, use
441 -pc_type ilu -pc_factor_mat_solver_type strumpack
442 to enable low-rank compression (i.e, use as a preconditioner).
443
444 Works with AIJ matrices
445
446 Options Database Keys:
447 + -mat_strumpack_verbose
448 . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
449 . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
450 . -mat_strumpack_colperm <TRUE> - Permute matrix to make diagonal nonzeros (None)
451 . -mat_strumpack_hss_min_sep_size <256> - Minimum size of separator for HSS compression (None)
452 . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
453 . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
454 . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
455 - -mat_strumpack_iterative_solver <DIRECT> - Select iterative solver from STRUMPACK (choose one of) AUTO DIRECT REFINE PREC_GMRES GMRES PREC_BICGSTAB BICGSTAB (None)
456
457 Level: beginner
458
459 .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
460 M*/
MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat * F)461 static PetscErrorCode MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat *F)
462 {
463 Mat B;
464 PetscErrorCode ierr;
465 PetscInt M=A->rmap->N,N=A->cmap->N;
466 PetscBool verb,flg,set;
467 PetscReal ctol;
468 PetscInt hssminsize,max_rank,leaf_size;
469 STRUMPACK_SparseSolver *S;
470 STRUMPACK_INTERFACE iface;
471 STRUMPACK_REORDERING_STRATEGY ndcurrent,ndvalue;
472 STRUMPACK_KRYLOV_SOLVER itcurrent,itsolver;
473 const STRUMPACK_PRECISION table[2][2][2] =
474 {{{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64},
475 {STRUMPACK_FLOAT_64, STRUMPACK_DOUBLE_64}},
476 {{STRUMPACK_FLOATCOMPLEX, STRUMPACK_DOUBLECOMPLEX},
477 {STRUMPACK_FLOAT, STRUMPACK_DOUBLE}}};
478 const STRUMPACK_PRECISION prec =
479 table[(sizeof(PetscInt)==8)?0:1]
480 [(PETSC_SCALAR==PETSC_COMPLEX)?0:1]
481 [(PETSC_REAL==PETSC_FLOAT)?0:1];
482 const char *const STRUMPACKNDTypes[] =
483 {"NATURAL","METIS","PARMETIS","SCOTCH","PTSCOTCH","RCM","STRUMPACKNDTypes","",0};
484 const char *const SolverTypes[] =
485 {"AUTO","NONE","REFINE","PREC_GMRES","GMRES","PREC_BICGSTAB","BICGSTAB","SolverTypes","",0};
486
487 PetscFunctionBegin;
488 /* Create the factorization matrix */
489 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
490 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);CHKERRQ(ierr);
491 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
492 ierr = MatSeqAIJSetPreallocation(B,0,NULL);
493 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr);
494 if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
495 B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
496 B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
497 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
498 B->ops->view = MatView_STRUMPACK;
499 B->ops->destroy = MatDestroy_STRUMPACK;
500 B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
501 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_aij_strumpack);CHKERRQ(ierr);
502 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetReordering_C",MatSTRUMPACKSetReordering_STRUMPACK);CHKERRQ(ierr);
503 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetColPerm_C",MatSTRUMPACKSetColPerm_STRUMPACK);CHKERRQ(ierr);
504 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSRelTol_C",MatSTRUMPACKSetHSSRelTol_STRUMPACK);CHKERRQ(ierr);
505 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSAbsTol_C",MatSTRUMPACKSetHSSAbsTol_STRUMPACK);CHKERRQ(ierr);
506 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMaxRank_C",MatSTRUMPACKSetHSSMaxRank_STRUMPACK);CHKERRQ(ierr);
507 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSLeafSize_C",MatSTRUMPACKSetHSSLeafSize_STRUMPACK);CHKERRQ(ierr);
508 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMinSepSize_C",MatSTRUMPACKSetHSSMinSepSize_STRUMPACK);CHKERRQ(ierr);
509 B->factortype = ftype;
510 ierr = PetscNewLog(B,&S);CHKERRQ(ierr);
511 B->spptr = S;
512
513 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
514 iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
515
516 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK Options","Mat");CHKERRQ(ierr);
517
518 verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
519 ierr = PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);CHKERRQ(ierr);
520
521 PetscStackCall("STRUMPACK_init",STRUMPACK_init(S,PetscObjectComm((PetscObject)A),prec,iface,0,NULL,verb));
522
523 PetscStackCall("STRUMPACK_HSS_rel_tol",ctol = (PetscReal)STRUMPACK_HSS_rel_tol(*S));
524 ierr = PetscOptionsReal("-mat_strumpack_hss_rel_tol","Relative compression tolerance","None",ctol,&ctol,&set);CHKERRQ(ierr);
525 if (set) PetscStackCall("STRUMPACK_set_HSS_rel_tol",STRUMPACK_set_HSS_rel_tol(*S,(double)ctol));
526
527 PetscStackCall("STRUMPACK_HSS_abs_tol",ctol = (PetscReal)STRUMPACK_HSS_abs_tol(*S));
528 ierr = PetscOptionsReal("-mat_strumpack_hss_abs_tol","Absolute compression tolerance","None",ctol,&ctol,&set);CHKERRQ(ierr);
529 if (set) PetscStackCall("STRUMPACK_set_HSS_abs_tol",STRUMPACK_set_HSS_abs_tol(*S,(double)ctol));
530
531 PetscStackCall("STRUMPACK_mc64job",flg = (STRUMPACK_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
532 ierr = PetscOptionsBool("-mat_strumpack_colperm","Find a col perm to get nonzero diagonal","None",flg,&flg,&set);CHKERRQ(ierr);
533 if (set) PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,flg ? 5 : 0));
534
535 PetscStackCall("STRUMPACK_HSS_min_sep_size",hssminsize = (PetscInt)STRUMPACK_HSS_min_sep_size(*S));
536 ierr = PetscOptionsInt("-mat_strumpack_hss_min_sep_size","Minimum size of separator for HSS compression","None",hssminsize,&hssminsize,&set);CHKERRQ(ierr);
537 if (set) PetscStackCall("STRUMPACK_set_HSS_min_sep_size",STRUMPACK_set_HSS_min_sep_size(*S,(int)hssminsize));
538
539 PetscStackCall("STRUMPACK_HSS_max_rank",max_rank = (PetscInt)STRUMPACK_HSS_max_rank(*S));
540 ierr = PetscOptionsInt("-mat_strumpack_max_rank","Maximum rank in HSS approximation","None",max_rank,&max_rank,&set);CHKERRQ(ierr);
541 if (set) PetscStackCall("STRUMPACK_set_HSS_max_rank",STRUMPACK_set_HSS_max_rank(*S,(int)max_rank));
542
543 PetscStackCall("STRUMPACK_HSS_leaf_size",leaf_size = (PetscInt)STRUMPACK_HSS_leaf_size(*S));
544 ierr = PetscOptionsInt("-mat_strumpack_leaf_size","Size of diagonal blocks in HSS approximation","None",leaf_size,&leaf_size,&set);CHKERRQ(ierr);
545 if (set) PetscStackCall("STRUMPACK_set_HSS_leaf_size",STRUMPACK_set_HSS_leaf_size(*S,(int)leaf_size));
546
547 PetscStackCall("STRUMPACK_reordering_method",ndcurrent = STRUMPACK_reordering_method(*S));
548 PetscOptionsEnum("-mat_strumpack_reordering","Sparsity reducing matrix reordering","None",STRUMPACKNDTypes,(PetscEnum)ndcurrent,(PetscEnum*)&ndvalue,&set);CHKERRQ(ierr);
549 if (set) PetscStackCall("STRUMPACK_set_reordering_method",STRUMPACK_set_reordering_method(*S,ndvalue));
550
551 if (ftype == MAT_FACTOR_ILU) {
552 /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete */
553 /* (or approximate) LU factorization. */
554 PetscStackCall("STRUMPACK_enable_HSS",STRUMPACK_enable_HSS(*S));
555 }
556
557 /* Disable the outer iterative solver from STRUMPACK. */
558 /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
559 /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling */
560 /* low-rank compression), it will use it's own preconditioned GMRES. Here we can disable */
561 /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
562 PetscStackCall("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
563
564 PetscStackCall("STRUMPACK_Krylov_solver",itcurrent = STRUMPACK_Krylov_solver(*S));
565 PetscOptionsEnum("-mat_strumpack_iterative_solver","Select iterative solver from STRUMPACK","None",SolverTypes,(PetscEnum)itcurrent,(PetscEnum*)&itsolver,&set);CHKERRQ(ierr);
566 if (set) PetscStackCall("STRUMPACK_set_Krylov_solver",STRUMPACK_set_Krylov_solver(*S,itsolver));
567
568 PetscOptionsEnd();
569
570 *F = B;
571 PetscFunctionReturn(0);
572 }
573
MatSolverTypeRegister_STRUMPACK(void)574 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
575 {
576 PetscErrorCode ierr;
577
578 PetscFunctionBegin;
579 ierr = MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
580 ierr = MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
581 ierr = MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
582 ierr = MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
583 PetscFunctionReturn(0);
584 }
585