xref: /petsc/src/mat/impls/aij/seq/umfpack/umfpack.c (revision 7a4fe282d1b349e95b3be72d69d8dd3d3bcd7bc6)
1 #define PETSCMAT_DLL
2 
3 /*
4    Provides an interface to the UMFPACKv5.1 sparse solver
5 
6    When build with PETSC_USE_64BIT_INDICES this will use UF_Long as the
7    integer type in UMFPACK, otherwise it will use int. This means
8    all integers in this file as simply declared as PetscInt. Also it means
9    that UMFPACK UL_Long version MUST be built with 64 bit integers when used.
10 
11 */
12 #include "../src/mat/impls/aij/seq/aij.h"
13 
14 #if defined(PETSC_USE_64BIT_INDICES)
15 #if defined(PETSC_USE_COMPLEX)
16 #define umfpack_UMF_free_symbolic   umfpack_zl_free_symbolic
17 #define umfpack_UMF_free_numeric    umfpack_zl_free_numeric
18 #define umfpack_UMF_wsolve          umfpack_zl_wsolve
19 #define umfpack_UMF_numeric         umfpack_zl_numeric
20 #define umfpack_UMF_report_numeric  umfpack_zl_report_numeric
21 #define umfpack_UMF_report_control  umfpack_zl_report_control
22 #define umfpack_UMF_report_status   umfpack_zl_report_status
23 #define umfpack_UMF_report_info     umfpack_zl_report_info
24 #define umfpack_UMF_report_symbolic umfpack_zl_report_symbolic
25 #define umfpack_UMF_qsymbolic       umfpack_zl_qsymbolic
26 #define umfpack_UMF_symbolic        umfpack_zl_symbolic
27 #define umfpack_UMF_defaults        umfpack_zl_defaults
28 
29 #else
30 #define umfpack_UMF_free_symbolic   umfpack_dl_free_symbolic
31 #define umfpack_UMF_free_numeric    umfpack_dl_free_numeric
32 #define umfpack_UMF_wsolve          umfpack_dl_wsolve
33 #define umfpack_UMF_numeric         umfpack_dl_numeric
34 #define umfpack_UMF_report_numeric  umfpack_dl_report_numeric
35 #define umfpack_UMF_report_control  umfpack_dl_report_control
36 #define umfpack_UMF_report_status   umfpack_dl_report_status
37 #define umfpack_UMF_report_info     umfpack_dl_report_info
38 #define umfpack_UMF_report_symbolic umfpack_dl_report_symbolic
39 #define umfpack_UMF_qsymbolic       umfpack_dl_qsymbolic
40 #define umfpack_UMF_symbolic        umfpack_dl_symbolic
41 #define umfpack_UMF_defaults        umfpack_dl_defaults
42 #endif
43 
44 #else
45 #if defined(PETSC_USE_COMPLEX)
46 #define umfpack_UMF_free_symbolic   umfpack_zi_free_symbolic
47 #define umfpack_UMF_free_numeric    umfpack_zi_free_numeric
48 #define umfpack_UMF_wsolve          umfpack_zi_wsolve
49 #define umfpack_UMF_numeric         umfpack_zi_numeric
50 #define umfpack_UMF_report_numeric  umfpack_zi_report_numeric
51 #define umfpack_UMF_report_control  umfpack_zi_report_control
52 #define umfpack_UMF_report_status   umfpack_zi_report_status
53 #define umfpack_UMF_report_info     umfpack_zi_report_info
54 #define umfpack_UMF_report_symbolic umfpack_zi_report_symbolic
55 #define umfpack_UMF_qsymbolic       umfpack_zi_qsymbolic
56 #define umfpack_UMF_symbolic        umfpack_zi_symbolic
57 #define umfpack_UMF_defaults        umfpack_zi_defaults
58 
59 #else
60 #define umfpack_UMF_free_symbolic   umfpack_di_free_symbolic
61 #define umfpack_UMF_free_numeric    umfpack_di_free_numeric
62 #define umfpack_UMF_wsolve          umfpack_di_wsolve
63 #define umfpack_UMF_numeric         umfpack_di_numeric
64 #define umfpack_UMF_report_numeric  umfpack_di_report_numeric
65 #define umfpack_UMF_report_control  umfpack_di_report_control
66 #define umfpack_UMF_report_status   umfpack_di_report_status
67 #define umfpack_UMF_report_info     umfpack_di_report_info
68 #define umfpack_UMF_report_symbolic umfpack_di_report_symbolic
69 #define umfpack_UMF_qsymbolic       umfpack_di_qsymbolic
70 #define umfpack_UMF_symbolic        umfpack_di_symbolic
71 #define umfpack_UMF_defaults        umfpack_di_defaults
72 #endif
73 #endif
74 
75 
76 #define UF_long long long
77 #define UF_long_max LONG_LONG_MAX
78 #define UF_long_id "%lld"
79 
80 EXTERN_C_BEGIN
81 #include "umfpack.h"
82 EXTERN_C_END
83 
84 typedef struct {
85   void         *Symbolic, *Numeric;
86   double       Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W;
87   PetscInt      *Wi,*ai,*aj,*perm_c;
88   PetscScalar  *av;
89   MatStructure flg;
90   PetscTruth   PetscMatOdering;
91 
92   /* Flag to clean up UMFPACK objects during Destroy */
93   PetscTruth CleanUpUMFPACK;
94 } Mat_UMFPACK;
95 
96 #undef __FUNCT__
97 #define __FUNCT__ "MatDestroy_UMFPACK"
98 PetscErrorCode MatDestroy_UMFPACK(Mat A)
99 {
100   PetscErrorCode ierr;
101   Mat_UMFPACK    *lu=(Mat_UMFPACK*)A->spptr;
102 
103   PetscFunctionBegin;
104   if (lu->CleanUpUMFPACK) {
105     umfpack_UMF_free_symbolic(&lu->Symbolic);
106     umfpack_UMF_free_numeric(&lu->Numeric);
107     ierr = PetscFree(lu->Wi);CHKERRQ(ierr);
108     ierr = PetscFree(lu->W);CHKERRQ(ierr);
109     if (lu->PetscMatOdering) {
110       ierr = PetscFree(lu->perm_c);CHKERRQ(ierr);
111     }
112   }
113   ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr);
114   PetscFunctionReturn(0);
115 }
116 
117 #undef __FUNCT__
118 #define __FUNCT__ "MatSolve_UMFPACK"
119 PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x)
120 {
121   Mat_UMFPACK    *lu = (Mat_UMFPACK*)A->spptr;
122   PetscScalar    *av=lu->av,*ba,*xa;
123   PetscErrorCode ierr;
124   PetscInt       *ai=lu->ai,*aj=lu->aj,status;
125 
126   PetscFunctionBegin;
127   /* solve Ax = b by umfpack_*_wsolve */
128   /* ----------------------------------*/
129   ierr = VecConjugate(b);
130 
131   ierr = VecGetArray(b,&ba);
132   ierr = VecGetArray(x,&xa);
133 #if defined(PETSC_USE_COMPLEX)
134   status = umfpack_UMF_wsolve(UMFPACK_At,ai,aj,(PetscReal*)av,NULL,(PetscReal*)xa,NULL,(PetscReal*)ba,NULL,
135                               lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
136 #else
137   status = umfpack_UMF_wsolve(UMFPACK_At,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
138 #endif
139   umfpack_UMF_report_info(lu->Control, lu->Info);
140   if (status < 0){
141     umfpack_UMF_report_status(lu->Control, status);
142     SETERRQ(PETSC_ERR_LIB,"umfpack_UMF_wsolve failed");
143   }
144 
145   ierr = VecRestoreArray(b,&ba);
146   ierr = VecRestoreArray(x,&xa);
147 
148   ierr = VecConjugate(b);
149   ierr = VecConjugate(x);
150   PetscFunctionReturn(0);
151 }
152 
153 #undef __FUNCT__
154 #define __FUNCT__ "MatLUFactorNumeric_UMFPACK"
155 PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info)
156 {
157   Mat_UMFPACK *lu=(Mat_UMFPACK*)(F)->spptr;
158   PetscErrorCode ierr;
159   PetscInt     *ai=lu->ai,*aj=lu->aj,m=A->rmap->n,status;
160   PetscScalar *av=lu->av;
161 
162   PetscFunctionBegin;
163   /* numeric factorization of A' */
164   /* ----------------------------*/
165 
166   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric){
167     umfpack_UMF_free_numeric(&lu->Numeric);
168   }
169 #if defined(PETSC_USE_COMPLEX)
170   status = umfpack_UMF_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
171 #else
172   status = umfpack_UMF_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
173 #endif
174   if (status < 0) {
175     umfpack_UMF_report_status(lu->Control, status);
176     SETERRQ(PETSC_ERR_LIB,"umfpack_UMF_numeric failed");
177   }
178   /* report numeric factorization of A' when Control[PRL] > 3 */
179   (void) umfpack_UMF_report_numeric(lu->Numeric, lu->Control);
180 
181   if (lu->flg == DIFFERENT_NONZERO_PATTERN){  /* first numeric factorization */
182     /* allocate working space to be used by Solve */
183     ierr = PetscMalloc(m * sizeof(PetscInt), &lu->Wi);CHKERRQ(ierr);
184     ierr = PetscMalloc(5*m * sizeof(PetscScalar), &lu->W);CHKERRQ(ierr);
185   }
186 
187   lu->flg = SAME_NONZERO_PATTERN;
188   lu->CleanUpUMFPACK = PETSC_TRUE;
189   (F)->ops->solve            = MatSolve_UMFPACK;
190   PetscFunctionReturn(0);
191 }
192 
193 /*
194    Note the r permutation is ignored
195 */
196 #undef __FUNCT__
197 #define __FUNCT__ "MatLUFactorSymbolic_UMFPACK"
198 PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
199 {
200   Mat_SeqAIJ     *mat=(Mat_SeqAIJ*)A->data;
201   Mat_UMFPACK    *lu = (Mat_UMFPACK*)(F->spptr);
202   PetscErrorCode ierr;
203   PetscInt       i,m=A->rmap->n,n=A->cmap->n;
204   const PetscInt *ra;
205   PetscInt        status;
206   PetscScalar    *av=mat->a;
207 
208   PetscFunctionBegin;
209   if (lu->PetscMatOdering) {
210     ierr = ISGetIndices(r,&ra);CHKERRQ(ierr);
211     ierr = PetscMalloc(m*sizeof(PetscInt),&lu->perm_c);CHKERRQ(ierr);
212     /* we cannot simply memcpy on 64 bit archs */
213     for(i = 0; i < m; i++) lu->perm_c[i] = ra[i];
214     ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr);
215   }
216 
217   lu->ai = mat->i;
218   lu->aj = mat->j;
219 
220   /* print the control parameters */
221   if(lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control);
222 
223   /* symbolic factorization of A' */
224   /* ---------------------------------------------------------------------- */
225   if (lu->PetscMatOdering) { /* use Petsc row ordering */
226 #if !defined(PETSC_USE_COMPLEX)
227     status = umfpack_UMF_qsymbolic(n,m,lu->ai,lu->aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
228 #else
229     status = umfpack_UMF_qsymbolic(n,m,lu->ai,lu->aj,NULL,NULL,
230                                    lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
231 #endif
232   } else { /* use Umfpack col ordering */
233 #if !defined(PETSC_USE_COMPLEX)
234     status = umfpack_UMF_symbolic(n,m,lu->ai,lu->aj,av,&lu->Symbolic,lu->Control,lu->Info);
235 #else
236     status = umfpack_UMF_symbolic(n,m,lu->ai,lu->aj,NULL,NULL,&lu->Symbolic,lu->Control,lu->Info);
237 #endif
238   }
239   if (status < 0){
240     umfpack_UMF_report_info(lu->Control, lu->Info);
241     umfpack_UMF_report_status(lu->Control, status);
242     SETERRQ(PETSC_ERR_LIB,"umfpack_UMF_symbolic failed");
243   }
244   /* report sumbolic factorization of A' when Control[PRL] > 3 */
245   (void) umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control);
246 
247   lu->flg = DIFFERENT_NONZERO_PATTERN;
248   lu->av  = av;
249   lu->CleanUpUMFPACK = PETSC_TRUE;
250   (F)->ops->lufactornumeric  = MatLUFactorNumeric_UMFPACK;
251   PetscFunctionReturn(0);
252 }
253 
254 #undef __FUNCT__
255 #define __FUNCT__ "MatFactorInfo_UMFPACK"
256 PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer)
257 {
258   Mat_UMFPACK    *lu= (Mat_UMFPACK*)A->spptr;
259   PetscErrorCode ierr;
260 
261   PetscFunctionBegin;
262   /* check if matrix is UMFPACK type */
263   if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0);
264 
265   ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr);
266   /* Control parameters used by reporting routiones */
267   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr);
268 
269   /* Control parameters used by symbolic factorization */
270   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);CHKERRQ(ierr);
271   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr);
272   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr);
273   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);CHKERRQ(ierr);
274   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr);
275   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_2BY2_TOLERANCE]: %g\n",lu->Control[UMFPACK_2BY2_TOLERANCE]);CHKERRQ(ierr);
276   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);CHKERRQ(ierr);
277   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);CHKERRQ(ierr);
278 
279   /* Control parameters used by numeric factorization */
280   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr);
281   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);CHKERRQ(ierr);
282   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);CHKERRQ(ierr);
283   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr);
284   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);CHKERRQ(ierr);
285 
286   /* Control parameters used by solve */
287   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr);
288 
289   /* mat ordering */
290   if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer,"  UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr);
291 
292   PetscFunctionReturn(0);
293 }
294 
295 #undef __FUNCT__
296 #define __FUNCT__ "MatView_UMFPACK"
297 PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer)
298 {
299   PetscErrorCode    ierr;
300   PetscTruth        iascii;
301   PetscViewerFormat format;
302 
303   PetscFunctionBegin;
304   ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr);
305 
306   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
307   if (iascii) {
308     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
309     if (format == PETSC_VIEWER_ASCII_INFO) {
310       ierr = MatFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr);
311     }
312   }
313   PetscFunctionReturn(0);
314 }
315 
316 EXTERN_C_BEGIN
317 #undef __FUNCT__
318 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_umfpack"
319 PetscErrorCode MatFactorGetSolverPackage_seqaij_umfpack(Mat A,const MatSolverPackage *type)
320 {
321   PetscFunctionBegin;
322   *type = MAT_SOLVER_UMFPACK;
323   PetscFunctionReturn(0);
324 }
325 EXTERN_C_END
326 
327 
328 /*MC
329   MAT_SOLVER_UMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices
330   via the external package UMFPACK.
331 
332   config/configure.py --download-umfpack to install PETSc to use UMFPACK
333 
334   Consult UMFPACK documentation for more information about the Control parameters
335   which correspond to the options database keys below.
336 
337   Options Database Keys:
338 + -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL]
339 . -mat_umfpack_strategy <AUTO> (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2
340 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
341 . -mat_umfpack_dense_row <0.2>: Control[UMFPACK_DENSE_ROW]
342 . -mat_umfpack_amd_dense <10>: Control[UMFPACK_AMD_DENSE]
343 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
344 . -mat_umfpack_2by2_tolerance <0.01>: Control[UMFPACK_2BY2_TOLERANCE]
345 . -mat_umfpack_fixq <0>: Control[UMFPACK_FIXQ]
346 . -mat_umfpack_aggressive <1>: Control[UMFPACK_AGGRESSIVE]
347 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
348 .  -mat_umfpack_sym_pivot_tolerance <0.001>: Control[UMFPACK_SYM_PIVOT_TOLERANCE]
349 .  -mat_umfpack_scale <NONE> (choose one of) NONE SUM MAX
350 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
351 .  -mat_umfpack_droptol <0>: Control[UMFPACK_DROPTOL]
352 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]
353 
354    Level: beginner
355 
356 .seealso: PCLU, MAT_SOLVER_SUPERLU, MAT_SOLVER_MUMPS, MAT_SOLVER_SPOOLES, PCFactorSetMatSolverPackage(), MatSolverPackage
357 M*/
358 EXTERN_C_BEGIN
359 #undef __FUNCT__
360 #define __FUNCT__ "MatGetFactor_seqaij_umfpack"
361 PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F)
362 {
363   Mat            B;
364   Mat_UMFPACK    *lu;
365   PetscErrorCode ierr;
366   PetscInt       m=A->rmap->n,n=A->cmap->n,idx;
367 
368   const char     *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC","2BY2"},
369                  *scale[]={"NONE","SUM","MAX"};
370   PetscTruth     flg;
371 
372   PetscFunctionBegin;
373   /* Create the factorization matrix F */
374   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
375   ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
376   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
377   ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
378   ierr = PetscNewLog(B,Mat_UMFPACK,&lu);CHKERRQ(ierr);
379   B->spptr                 = lu;
380   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
381   B->ops->destroy          = MatDestroy_UMFPACK;
382   B->ops->view             = MatView_UMFPACK;
383   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_umfpack",MatFactorGetSolverPackage_seqaij_umfpack);CHKERRQ(ierr);
384   B->factor                = MAT_FACTOR_LU;
385   B->assembled             = PETSC_TRUE;  /* required by -ksp_view */
386   B->preallocated          = PETSC_TRUE;
387 
388   /* initializations */
389   /* ------------------------------------------------*/
390   /* get the default control parameters */
391   umfpack_UMF_defaults(lu->Control);
392   lu->perm_c = PETSC_NULL;  /* use defaul UMFPACK col permutation */
393   lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */
394 
395   ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr);
396   /* Control parameters used by reporting routiones */
397   ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr);
398 
399   /* Control parameters for symbolic factorization */
400   ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,4,strategy[0],&idx,&flg);CHKERRQ(ierr);
401   if (flg) {
402     switch (idx){
403     case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break;
404     case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break;
405     case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break;
406     case 3: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2; break;
407     }
408   }
409   ierr = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],PETSC_NULL);CHKERRQ(ierr);
410   ierr = PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],PETSC_NULL);CHKERRQ(ierr);
411   ierr = PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],PETSC_NULL);CHKERRQ(ierr);
412   ierr = PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],PETSC_NULL);CHKERRQ(ierr);
413   ierr = PetscOptionsReal("-mat_umfpack_2by2_tolerance","Control[UMFPACK_2BY2_TOLERANCE]","None",lu->Control[UMFPACK_2BY2_TOLERANCE],&lu->Control[UMFPACK_2BY2_TOLERANCE],PETSC_NULL);CHKERRQ(ierr);
414   ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);CHKERRQ(ierr);
415   ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);CHKERRQ(ierr);
416 
417   /* Control parameters used by numeric factorization */
418   ierr = PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],PETSC_NULL);CHKERRQ(ierr);
419   ierr = PetscOptionsReal("-mat_umfpack_sym_pivot_tolerance","Control[UMFPACK_SYM_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],&lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],PETSC_NULL);CHKERRQ(ierr);
420   ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr);
421   if (flg) {
422     switch (idx){
423     case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
424     case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
425     case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break;
426     }
427   }
428   ierr = PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);CHKERRQ(ierr);
429   ierr = PetscOptionsReal("-mat_umfpack_front_alloc_init","Control[UMFPACK_FRONT_ALLOC_INIT]","None",lu->Control[UMFPACK_FRONT_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);CHKERRQ(ierr);
430   ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],PETSC_NULL);CHKERRQ(ierr);
431 
432   /* Control parameters used by solve */
433   ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr);
434 
435   /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */
436   ierr = PetscOptionsHasName(PETSC_NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr);
437   PetscOptionsEnd();
438   *F = B;
439   PetscFunctionReturn(0);
440 }
441 EXTERN_C_END
442 
443