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