xref: /petsc/src/mat/impls/aij/mpi/strumpack/strumpack.c (revision 4fc747eaadbeca11629f314a99edccbc2ed7b3d3)
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 
5 #undef __FUNCT__
6 #define __FUNCT__ "MatGetDiagonal_STRUMPACK"
7 static PetscErrorCode MatGetDiagonal_STRUMPACK(Mat A,Vec v)
8 {
9   PetscFunctionBegin;
10   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK factor");
11   PetscFunctionReturn(0);
12 }
13 
14 #undef __FUNCT__
15 #define __FUNCT__ "MatDestroy_STRUMPACK"
16 static PetscErrorCode MatDestroy_STRUMPACK(Mat A)
17 {
18   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
19   PetscErrorCode         ierr;
20   PetscBool              flg;
21 
22   PetscFunctionBegin;
23   /* Deallocate STRUMPACK storage */
24   PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(S));
25   ierr = PetscFree(A->spptr);CHKERRQ(ierr);
26   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr);
27   if (flg) {
28     ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr);
29   } else {
30     ierr = MatDestroy_MPIAIJ(A);CHKERRQ(ierr);
31   }
32 
33   /* clear composed functions */
34   ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr);
35   ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSRelCompTol_C",NULL);CHKERRQ(ierr);
36   ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMinSize_C",NULL);CHKERRQ(ierr);
37   ierr = PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetColPerm_C",NULL);CHKERRQ(ierr);
38 
39   PetscFunctionReturn(0);
40 }
41 
42 #undef __FUNCT__
43 #define __FUNCT__ "MatSTRUMPACKSetHSSRelCompTol_STRUMPACK"
44 static PetscErrorCode MatSTRUMPACKSetHSSRelCompTol_STRUMPACK(Mat F,PetscReal rctol)
45 {
46   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
47 
48   PetscFunctionBegin;
49   PetscStackCall("STRUMPACK_set_rctol", STRUMPACK_set_rctol(*S,rctol));
50   PetscFunctionReturn(0);
51 }
52 
53 #undef __FUNCT__
54 #define __FUNCT__ "MatSTRUMPACKSetHSSRelCompTol"
55 /*@
56   MatSTRUMPACKSetHSSRelCompTol - Set STRUMPACK relative tolerance for HSS compression
57    Logically Collective on Mat
58 
59    Input Parameters:
60 +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
61 -  rctol - relative compression tolerance
62 
63   Options Database:
64 .   -mat_strumpack_rctol <rctol>
65 
66    Level: beginner
67 
68    References:
69 .      STRUMPACK manual
70 
71 .seealso: MatGetFactor()
72 @*/
73 PetscErrorCode MatSTRUMPACKSetHSSRelCompTol(Mat F,PetscReal rctol)
74 {
75   PetscErrorCode ierr;
76 
77   PetscFunctionBegin;
78   PetscValidHeaderSpecific(F,MAT_CLASSID,1);
79   PetscValidLogicalCollectiveReal(F,rctol,2);
80   ierr = PetscTryMethod(F,"MatSTRUMPACKSetHSSRelCompTol_C",(Mat,PetscReal),(F,rctol));CHKERRQ(ierr);
81   PetscFunctionReturn(0);
82 }
83 
84 #undef __FUNCT__
85 #define __FUNCT__ "MatSTRUMPACKSetHSSMinSize_STRUMPACK"
86 static PetscErrorCode MatSTRUMPACKSetHSSMinSize_STRUMPACK(Mat F,PetscInt hssminsize)
87 {
88   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
89 
90   PetscFunctionBegin;
91   PetscStackCall("STRUMPACK_set_minimum_HSS_size", STRUMPACK_set_minimum_HSS_size(*S,hssminsize));
92   PetscFunctionReturn(0);
93 }
94 
95 #undef __FUNCT__
96 #define __FUNCT__ "MatSTRUMPACKSetHSSMinSize"
97 /*@
98   MatSTRUMPACKSetHSSMinSize - Set STRUMPACK minimum dense matrix size for low-rank approximation
99    Logically Collective on Mat
100 
101    Input Parameters:
102 +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
103 -  hssminsize - minimum dense matrix size for low-rank approximation
104 
105   Options Database:
106 .   -mat_strumpack_hssminsize <hssminsize>
107 
108    Level: beginner
109 
110    References:
111 .      STRUMPACK manual
112 
113 .seealso: MatGetFactor()
114 @*/
115 PetscErrorCode MatSTRUMPACKSetHSSMinSize(Mat F,PetscInt hssminsize)
116 {
117   PetscErrorCode ierr;
118 
119   PetscFunctionBegin;
120   PetscValidHeaderSpecific(F,MAT_CLASSID,1);
121   PetscValidLogicalCollectiveInt(F,hssminsize,2);
122   ierr = PetscTryMethod(F,"MatSTRUMPACKSetHSSMinSize_C",(Mat,PetscInt),(F,hssminsize));CHKERRQ(ierr);
123   PetscFunctionReturn(0);
124 }
125 
126 #undef __FUNCT__
127 #define __FUNCT__ "MatSTRUMPACKSetColPerm_STRUMPACK"
128 static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)
129 {
130   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
131 
132   PetscFunctionBegin;
133   PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,cperm ? 5 : 0));
134   PetscFunctionReturn(0);
135 }
136 
137 #undef __FUNCT__
138 #define __FUNCT__ "MatSTRUMPACKSetColPerm"
139 /*@
140   MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
141    Logically Collective on Mat
142 
143    Input Parameters:
144 +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
145 -  cperm - whether or not to permute (internally) the columns of the matrix
146 
147   Options Database:
148 .   -mat_strumpack_colperm <cperm>
149 
150    Level: beginner
151 
152    References:
153 .      STRUMPACK manual
154 
155 .seealso: MatGetFactor()
156 @*/
157 PetscErrorCode MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)
158 {
159   PetscErrorCode ierr;
160 
161   PetscFunctionBegin;
162   PetscValidHeaderSpecific(F,MAT_CLASSID,1);
163   PetscValidLogicalCollectiveBool(F,cperm,2);
164   ierr = PetscTryMethod(F,"MatSTRUMPACKSetColPerm_C",(Mat,PetscBool),(F,cperm));CHKERRQ(ierr);
165   PetscFunctionReturn(0);
166 }
167 
168 #undef __FUNCT__
169 #define __FUNCT__ "MatSolve_STRUMPACK"
170 static PetscErrorCode MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)
171 {
172   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
173   STRUMPACK_RETURN_CODE  sp_err;
174   PetscErrorCode         ierr;
175   const PetscScalar      *bptr;
176   PetscScalar            *xptr;
177 
178   PetscFunctionBegin;
179   ierr = VecGetArray(x,&xptr);CHKERRQ(ierr);
180   ierr = VecGetArrayRead(b_mpi,&bptr);CHKERRQ(ierr);
181 
182   PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(*S,(PetscScalar*)bptr,xptr,0));
183   switch (sp_err) {
184   case STRUMPACK_SUCCESS: break;
185   case STRUMPACK_MATRIX_NOT_SET:   { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
186   case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
187   default:                           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
188   }
189   ierr = VecRestoreArray(x,&xptr);CHKERRQ(ierr);
190   ierr = VecRestoreArrayRead(b_mpi,&bptr);CHKERRQ(ierr);
191   PetscFunctionReturn(0);
192 }
193 
194 #undef __FUNCT__
195 #define __FUNCT__ "MatMatSolve_STRUMPACK"
196 static PetscErrorCode MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)
197 {
198   PetscErrorCode   ierr;
199   PetscBool        flg;
200 
201   PetscFunctionBegin;
202   ierr = PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr);
203   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
204   ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr);
205   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
206   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK() is not implemented yet");
207   PetscFunctionReturn(0);
208 }
209 
210 #undef __FUNCT__
211 #define __FUNCT__ "MatFactorInfo_STRUMPACK"
212 static PetscErrorCode MatFactorInfo_STRUMPACK(Mat A,PetscViewer viewer)
213 {
214   PetscErrorCode  ierr;
215 
216   PetscFunctionBegin;
217   /* check if matrix is strumpack type */
218   if (A->ops->solve != MatSolve_STRUMPACK) PetscFunctionReturn(0);
219   ierr = PetscViewerASCIIPrintf(viewer,"STRUMPACK sparse solver!\n");CHKERRQ(ierr);
220   PetscFunctionReturn(0);
221 }
222 
223 #undef __FUNCT__
224 #define __FUNCT__ "MatView_STRUMPACK"
225 static PetscErrorCode MatView_STRUMPACK(Mat A,PetscViewer viewer)
226 {
227   PetscErrorCode    ierr;
228   PetscBool         iascii;
229   PetscViewerFormat format;
230 
231   PetscFunctionBegin;
232   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
233   if (iascii) {
234     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
235     if (format == PETSC_VIEWER_ASCII_INFO) {
236       ierr = MatFactorInfo_STRUMPACK(A,viewer);CHKERRQ(ierr);
237     }
238   }
239   PetscFunctionReturn(0);
240 }
241 
242 #undef __FUNCT__
243 #define __FUNCT__ "MatLUFactorNumeric_STRUMPACK"
244 static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo *info)
245 {
246   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
247   STRUMPACK_RETURN_CODE  sp_err;
248   Mat_SeqAIJ             *A_d,*A_o;
249   Mat_MPIAIJ             *mat;
250   PetscErrorCode         ierr;
251   PetscInt               M=A->rmap->N,m=A->rmap->n;
252   PetscBool              flg;
253 
254   PetscFunctionBegin;
255   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr);
256   if (flg) { /* A is MATMPIAIJ */
257     mat = (Mat_MPIAIJ*)A->data;
258     A_d = (Mat_SeqAIJ*)(mat->A)->data;
259     A_o = (Mat_SeqAIJ*)(mat->B)->data;
260     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));
261   } else { /* A is MATSEQAIJ */
262     A_d = (Mat_SeqAIJ*)A->data;
263     PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(*S,&M,A_d->i,A_d->j,A_d->a,0));
264   }
265 
266   /* Reorder and Factor the matrix. */
267   /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
268   PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(*S));
269   PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(*S));
270   switch (sp_err) {
271   case STRUMPACK_SUCCESS: break;
272   case STRUMPACK_MATRIX_NOT_SET:   { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
273   case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
274   default:                           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
275   }
276   PetscFunctionReturn(0);
277 }
278 
279 #undef __FUNCT__
280 #define __FUNCT__ "MatLUFactorSymbolic_STRUMPACK"
281 static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
282 {
283   PetscFunctionBegin;
284   F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
285   F->ops->solve           = MatSolve_STRUMPACK;
286   F->ops->matsolve        = MatMatSolve_STRUMPACK;
287   PetscFunctionReturn(0);
288 }
289 
290 #undef __FUNCT__
291 #define __FUNCT__ "MatFactorGetSolverPackage_aij_strumpack"
292 static PetscErrorCode MatFactorGetSolverPackage_aij_strumpack(Mat A,const MatSolverPackage *type)
293 {
294   PetscFunctionBegin;
295   *type = MATSOLVERSTRUMPACK;
296   PetscFunctionReturn(0);
297 }
298 
299 /*MC
300   MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
301   and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.
302 
303   Consult the STRUMPACK-sparse manual for more info.
304 
305   Use
306      ./configure --download-strumpack
307   to have PETSc installed with STRUMPACK
308 
309   Use
310     -pc_type lu -pc_factor_mat_solver_package strumpack
311   to use this as an exact (direct) solver, use
312     -pc_type ilu -pc_factor_mat_solver_package strumpack
313   to enable low-rank compression (i.e, use as a preconditioner).
314 
315   Works with AIJ matrices
316 
317   Options Database Keys:
318 + -mat_strumpack_rctol <1e-4>           - Relative compression tolerance (None)
319 . -mat_strumpack_hssminsize <512>       - Minimum size of dense block for HSS compression (None)
320 . -mat_strumpack_colperm <TRUE>         - Permute matrix to make diagonal nonzeros (None)
321 
322  Level: beginner
323 
324 .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverPackage(), MatSolverPackage
325 M*/
326 #undef __FUNCT__
327 #define __FUNCT__ "MatGetFactor_aij_strumpack"
328 static PetscErrorCode MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat *F)
329 {
330   Mat                    B;
331   STRUMPACK_SparseSolver *S;
332   PetscErrorCode         ierr;
333   PetscInt               M=A->rmap->N,N=A->cmap->N;
334   STRUMPACK_INTERFACE    iface;
335   PetscBool              verb,flg,set;
336   PetscReal              rctol;
337   PetscInt               hssminsize;
338   int                    argc;
339   char                   **args,*copts,*pname;
340   size_t                 len;
341   const STRUMPACK_PRECISION table[2][2][2] = {{{STRUMPACK_FLOATCOMPLEX_64,STRUMPACK_DOUBLECOMPLEX_64},
342                                                {STRUMPACK_FLOAT_64,STRUMPACK_DOUBLE_64}},
343                                               {{STRUMPACK_FLOATCOMPLEX,STRUMPACK_DOUBLECOMPLEX},
344                                                {STRUMPACK_FLOAT,STRUMPACK_DOUBLE}}};
345   const STRUMPACK_PRECISION prec = table[(sizeof(PetscInt)==8)?0:1][(PETSC_SCALAR==PETSC_COMPLEX)?0:1][(PETSC_REAL==PETSC_FLOAT)?0:1];
346 
347   PetscFunctionBegin;
348   /* Create the factorization matrix */
349   ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
350   ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);CHKERRQ(ierr);
351   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
352   ierr = MatSeqAIJSetPreallocation(B,0,NULL);
353   ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr);
354   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
355     B->ops->lufactorsymbolic  = MatLUFactorSymbolic_STRUMPACK;
356     B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
357   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
358   B->ops->view        = MatView_STRUMPACK;
359   B->ops->destroy     = MatDestroy_STRUMPACK;
360   B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
361   ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_aij_strumpack);CHKERRQ(ierr);
362   ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSRelCompTol_C",MatSTRUMPACKSetHSSRelCompTol_STRUMPACK);CHKERRQ(ierr);
363   ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMinSize_C",MatSTRUMPACKSetHSSMinSize_STRUMPACK);CHKERRQ(ierr);
364   ierr = PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetColPerm_C",MatSTRUMPACKSetColPerm_STRUMPACK);CHKERRQ(ierr);
365   B->factortype = ftype;
366   ierr     = PetscNewLog(B,&S);CHKERRQ(ierr);
367   B->spptr = S;
368 
369   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
370   iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
371 
372   ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK Options","Mat");CHKERRQ(ierr);
373 
374   verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
375   ierr = PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);CHKERRQ(ierr);
376 
377   ierr = PetscOptionsGetAll(NULL,&copts);CHKERRQ(ierr);
378   ierr = PetscStrlen(copts,&len);CHKERRQ(ierr);
379   len += PETSC_MAX_PATH_LEN+1;
380   ierr = PetscMalloc1(len,&pname);CHKERRQ(ierr);
381   /* first string is assumed to be the program name, so add program name to options string */
382   ierr = PetscGetProgramName(pname,len);CHKERRQ(ierr);
383   ierr = PetscStrcat(pname," ");CHKERRQ(ierr);
384   ierr = PetscStrcat(pname,copts);CHKERRQ(ierr);
385   ierr = PetscStrToArray(pname,' ',&argc,&args);CHKERRQ(ierr);
386   ierr = PetscFree(copts);CHKERRQ(ierr);
387   ierr = PetscFree(pname);CHKERRQ(ierr);
388 
389   PetscStackCall("STRUMPACK_init",STRUMPACK_init(S,PetscObjectComm((PetscObject)A),prec,iface,argc,args,verb));
390 
391   PetscStackCall("STRUMPACK_get_rctol",rctol = (PetscReal)STRUMPACK_get_rctol(*S));
392   ierr = PetscOptionsReal("-mat_strumpack_rctol","Relative compression tolerance","None",rctol,&rctol,&set);CHKERRQ(ierr);
393   if (set) PetscStackCall("STRUMPACK_set_rctol",STRUMPACK_set_rctol(*S,(double)rctol));
394 
395   PetscStackCall("STRUMPACK_get_mc64job",flg = (STRUMPACK_get_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
396   ierr = PetscOptionsBool("-mat_strumpack_colperm","Find a col perm to get nonzero diagonal","None",flg,&flg,&set);CHKERRQ(ierr);
397   if (set) PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,flg ? 5 : 0));
398 
399   PetscStackCall("STRUMPACK_get_minimum_HSS_size",hssminsize = (PetscInt)STRUMPACK_get_minimum_HSS_size(*S));
400   ierr = PetscOptionsInt("-mat_strumpack_hssminsize","Minimum size of dense block for HSS compression","None",hssminsize,&hssminsize,&set);CHKERRQ(ierr);
401   if (set) PetscStackCall("STRUMPACK_set_minimum_HSS_size",STRUMPACK_set_minimum_HSS_size(*S,(int)hssminsize));
402 
403   PetscOptionsEnd();
404 
405   if (ftype == MAT_FACTOR_ILU) {
406     /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete                */
407     /* (or approximate) LU factorization.                                                       */
408     PetscStackCall("STRUMPACK_use_HSS",STRUMPACK_use_HSS(*S,1));
409     /* Disable the outer iterative solver from STRUMPACK.                                       */
410     /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement.   */
411     /* When STRUMPACK is used with as an approximate factorization preconditioner (by enabling  */
412     /* low-rank compression), it will use it's own GMRES. Here we can disable the               */
413     /* outer iterative solver, as PETSc uses STRUMPACK from within a KSP.                       */
414     PetscStackCall("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
415   }
416 
417   /* Allow the user to set or overwrite the above options from the command line                 */
418   PetscStackCall("STRUMPACK_set_from_options",STRUMPACK_set_from_options(*S));
419   ierr = PetscStrToArrayDestroy(argc,args);CHKERRQ(ierr);
420 
421   *F = B;
422   PetscFunctionReturn(0);
423 }
424 
425 #undef __FUNCT__
426 #define __FUNCT__ "MatSolverPackageRegister_STRUMPACK"
427 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_STRUMPACK(void)
428 {
429   PetscErrorCode ierr;
430 
431   PetscFunctionBegin;
432   ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
433   ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
434   ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
435   ierr = MatSolverPackageRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);CHKERRQ(ierr);
436   PetscFunctionReturn(0);
437 }
438