xref: /petsc/src/mat/impls/aij/seq/umfpack/umfpack.c (revision 0700a8246d308f50502909ba325e6169d3ee27eb)
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 static 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_Private"
119 static PetscErrorCode MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag)
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 
130   ierr = VecGetArray(b,&ba);
131   ierr = VecGetArray(x,&xa);
132 #if defined(PETSC_USE_COMPLEX)
133   status = umfpack_UMF_wsolve(uflag,ai,aj,(PetscReal*)av,NULL,(PetscReal*)xa,NULL,(PetscReal*)ba,NULL,
134                               lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
135 #else
136   status = umfpack_UMF_wsolve(uflag,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
137 #endif
138   umfpack_UMF_report_info(lu->Control, lu->Info);
139   if (status < 0){
140     umfpack_UMF_report_status(lu->Control, status);
141     SETERRQ(PETSC_ERR_LIB,"umfpack_UMF_wsolve failed");
142   }
143 
144   ierr = VecRestoreArray(b,&ba);
145   ierr = VecRestoreArray(x,&xa);
146   PetscFunctionReturn(0);
147 }
148 
149 #undef __FUNCT__
150 #define __FUNCT__ "MatSolve_UMFPACK"
151 static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x)
152 {
153   PetscErrorCode ierr;
154 
155   PetscFunctionBegin;
156   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
157   ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_Aat);CHKERRQ(ierr);
158   PetscFunctionReturn(0);
159 }
160 
161 #undef __FUNCT__
162 #define __FUNCT__ "MatSolveTranspose_UMFPACK"
163 static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x)
164 {
165   PetscErrorCode ierr;
166 
167   PetscFunctionBegin;
168   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
169   ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A);CHKERRQ(ierr);
170   PetscFunctionReturn(0);
171 }
172 
173 #undef __FUNCT__
174 #define __FUNCT__ "MatLUFactorNumeric_UMFPACK"
175 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info)
176 {
177   Mat_UMFPACK *lu=(Mat_UMFPACK*)(F)->spptr;
178   PetscErrorCode ierr;
179   PetscInt     *ai=lu->ai,*aj=lu->aj,m=A->rmap->n,status;
180   PetscScalar *av=lu->av;
181 
182   PetscFunctionBegin;
183   /* numeric factorization of A' */
184   /* ----------------------------*/
185 
186   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric){
187     umfpack_UMF_free_numeric(&lu->Numeric);
188   }
189 #if defined(PETSC_USE_COMPLEX)
190   status = umfpack_UMF_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
191 #else
192   status = umfpack_UMF_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
193 #endif
194   if (status < 0) {
195     umfpack_UMF_report_status(lu->Control, status);
196     SETERRQ(PETSC_ERR_LIB,"umfpack_UMF_numeric failed");
197   }
198   /* report numeric factorization of A' when Control[PRL] > 3 */
199   (void) umfpack_UMF_report_numeric(lu->Numeric, lu->Control);
200 
201   if (lu->flg == DIFFERENT_NONZERO_PATTERN){  /* first numeric factorization */
202     /* allocate working space to be used by Solve */
203     ierr = PetscMalloc(m * sizeof(PetscInt), &lu->Wi);CHKERRQ(ierr);
204     ierr = PetscMalloc(5*m * sizeof(PetscScalar), &lu->W);CHKERRQ(ierr);
205   }
206 
207   lu->flg = SAME_NONZERO_PATTERN;
208   lu->CleanUpUMFPACK = PETSC_TRUE;
209   F->ops->solve          = MatSolve_UMFPACK;
210   F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
211   PetscFunctionReturn(0);
212 }
213 
214 /*
215    Note the r permutation is ignored
216 */
217 #undef __FUNCT__
218 #define __FUNCT__ "MatLUFactorSymbolic_UMFPACK"
219 static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
220 {
221   Mat_SeqAIJ     *mat=(Mat_SeqAIJ*)A->data;
222   Mat_UMFPACK    *lu = (Mat_UMFPACK*)(F->spptr);
223   PetscErrorCode ierr;
224   PetscInt       i,m=A->rmap->n,n=A->cmap->n;
225   const PetscInt *ra;
226   PetscInt        status;
227   PetscScalar    *av=mat->a;
228 
229   PetscFunctionBegin;
230   if (lu->PetscMatOdering) {
231     ierr = ISGetIndices(r,&ra);CHKERRQ(ierr);
232     ierr = PetscMalloc(m*sizeof(PetscInt),&lu->perm_c);CHKERRQ(ierr);
233     /* we cannot simply memcpy on 64 bit archs */
234     for(i = 0; i < m; i++) lu->perm_c[i] = ra[i];
235     ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr);
236   }
237 
238   lu->ai = mat->i;
239   lu->aj = mat->j;
240 
241   /* print the control parameters */
242   if(lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control);
243 
244   /* symbolic factorization of A' */
245   /* ---------------------------------------------------------------------- */
246   if (lu->PetscMatOdering) { /* use Petsc row ordering */
247 #if !defined(PETSC_USE_COMPLEX)
248     status = umfpack_UMF_qsymbolic(n,m,lu->ai,lu->aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
249 #else
250     status = umfpack_UMF_qsymbolic(n,m,lu->ai,lu->aj,NULL,NULL,
251                                    lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
252 #endif
253   } else { /* use Umfpack col ordering */
254 #if !defined(PETSC_USE_COMPLEX)
255     status = umfpack_UMF_symbolic(n,m,lu->ai,lu->aj,av,&lu->Symbolic,lu->Control,lu->Info);
256 #else
257     status = umfpack_UMF_symbolic(n,m,lu->ai,lu->aj,NULL,NULL,&lu->Symbolic,lu->Control,lu->Info);
258 #endif
259   }
260   if (status < 0){
261     umfpack_UMF_report_info(lu->Control, lu->Info);
262     umfpack_UMF_report_status(lu->Control, status);
263     SETERRQ(PETSC_ERR_LIB,"umfpack_UMF_symbolic failed");
264   }
265   /* report sumbolic factorization of A' when Control[PRL] > 3 */
266   (void) umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control);
267 
268   lu->flg = DIFFERENT_NONZERO_PATTERN;
269   lu->av  = av;
270   lu->CleanUpUMFPACK = PETSC_TRUE;
271   (F)->ops->lufactornumeric  = MatLUFactorNumeric_UMFPACK;
272   PetscFunctionReturn(0);
273 }
274 
275 #undef __FUNCT__
276 #define __FUNCT__ "MatFactorInfo_UMFPACK"
277 static PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer)
278 {
279   Mat_UMFPACK    *lu= (Mat_UMFPACK*)A->spptr;
280   PetscErrorCode ierr;
281 
282   PetscFunctionBegin;
283   /* check if matrix is UMFPACK type */
284   if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0);
285 
286   ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr);
287   /* Control parameters used by reporting routiones */
288   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr);
289 
290   /* Control parameters used by symbolic factorization */
291   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);CHKERRQ(ierr);
292   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr);
293   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr);
294   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);CHKERRQ(ierr);
295   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr);
296   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_2BY2_TOLERANCE]: %g\n",lu->Control[UMFPACK_2BY2_TOLERANCE]);CHKERRQ(ierr);
297   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);CHKERRQ(ierr);
298   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);CHKERRQ(ierr);
299 
300   /* Control parameters used by numeric factorization */
301   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr);
302   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);CHKERRQ(ierr);
303   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);CHKERRQ(ierr);
304   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr);
305   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);CHKERRQ(ierr);
306 
307   /* Control parameters used by solve */
308   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr);
309 
310   /* mat ordering */
311   if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer,"  UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr);
312 
313   PetscFunctionReturn(0);
314 }
315 
316 #undef __FUNCT__
317 #define __FUNCT__ "MatView_UMFPACK"
318 static PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer)
319 {
320   PetscErrorCode    ierr;
321   PetscTruth        iascii;
322   PetscViewerFormat format;
323 
324   PetscFunctionBegin;
325   ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr);
326 
327   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
328   if (iascii) {
329     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
330     if (format == PETSC_VIEWER_ASCII_INFO) {
331       ierr = MatFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr);
332     }
333   }
334   PetscFunctionReturn(0);
335 }
336 
337 EXTERN_C_BEGIN
338 #undef __FUNCT__
339 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_umfpack"
340 PetscErrorCode MatFactorGetSolverPackage_seqaij_umfpack(Mat A,const MatSolverPackage *type)
341 {
342   PetscFunctionBegin;
343   *type = MAT_SOLVER_UMFPACK;
344   PetscFunctionReturn(0);
345 }
346 EXTERN_C_END
347 
348 
349 /*MC
350   MAT_SOLVER_UMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices
351   via the external package UMFPACK.
352 
353   config/configure.py --download-umfpack to install PETSc to use UMFPACK
354 
355   Consult UMFPACK documentation for more information about the Control parameters
356   which correspond to the options database keys below.
357 
358   Options Database Keys:
359 + -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL]
360 . -mat_umfpack_strategy <AUTO> (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2
361 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
362 . -mat_umfpack_dense_row <0.2>: Control[UMFPACK_DENSE_ROW]
363 . -mat_umfpack_amd_dense <10>: Control[UMFPACK_AMD_DENSE]
364 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
365 . -mat_umfpack_2by2_tolerance <0.01>: Control[UMFPACK_2BY2_TOLERANCE]
366 . -mat_umfpack_fixq <0>: Control[UMFPACK_FIXQ]
367 . -mat_umfpack_aggressive <1>: Control[UMFPACK_AGGRESSIVE]
368 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
369 .  -mat_umfpack_sym_pivot_tolerance <0.001>: Control[UMFPACK_SYM_PIVOT_TOLERANCE]
370 .  -mat_umfpack_scale <NONE> (choose one of) NONE SUM MAX
371 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
372 .  -mat_umfpack_droptol <0>: Control[UMFPACK_DROPTOL]
373 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]
374 
375    Level: beginner
376 
377 .seealso: PCLU, MAT_SOLVER_SUPERLU, MAT_SOLVER_MUMPS, MAT_SOLVER_SPOOLES, PCFactorSetMatSolverPackage(), MatSolverPackage
378 M*/
379 EXTERN_C_BEGIN
380 #undef __FUNCT__
381 #define __FUNCT__ "MatGetFactor_seqaij_umfpack"
382 PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F)
383 {
384   Mat            B;
385   Mat_UMFPACK    *lu;
386   PetscErrorCode ierr;
387   PetscInt       m=A->rmap->n,n=A->cmap->n,idx;
388 
389   const char     *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC","2BY2"},
390                  *scale[]={"NONE","SUM","MAX"};
391   PetscTruth     flg;
392 
393   PetscFunctionBegin;
394   /* Create the factorization matrix F */
395   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
396   ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
397   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
398   ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
399   ierr = PetscNewLog(B,Mat_UMFPACK,&lu);CHKERRQ(ierr);
400   B->spptr                 = lu;
401   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
402   B->ops->destroy          = MatDestroy_UMFPACK;
403   B->ops->view             = MatView_UMFPACK;
404   ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_umfpack",MatFactorGetSolverPackage_seqaij_umfpack);CHKERRQ(ierr);
405   B->factortype            = MAT_FACTOR_LU;
406   B->assembled             = PETSC_TRUE;  /* required by -ksp_view */
407   B->preallocated          = PETSC_TRUE;
408 
409   /* initializations */
410   /* ------------------------------------------------*/
411   /* get the default control parameters */
412   umfpack_UMF_defaults(lu->Control);
413   lu->perm_c = PETSC_NULL;  /* use defaul UMFPACK col permutation */
414   lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */
415 
416   ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr);
417   /* Control parameters used by reporting routiones */
418   ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr);
419 
420   /* Control parameters for symbolic factorization */
421   ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,4,strategy[0],&idx,&flg);CHKERRQ(ierr);
422   if (flg) {
423     switch (idx){
424     case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break;
425     case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break;
426     case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break;
427     case 3: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2; break;
428     }
429   }
430   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);
431   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);
432   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);
433   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);
434   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);
435   ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);CHKERRQ(ierr);
436   ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);CHKERRQ(ierr);
437 
438   /* Control parameters used by numeric factorization */
439   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);
440   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);
441   ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr);
442   if (flg) {
443     switch (idx){
444     case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
445     case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
446     case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break;
447     }
448   }
449   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);
450   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);
451   ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],PETSC_NULL);CHKERRQ(ierr);
452 
453   /* Control parameters used by solve */
454   ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr);
455 
456   /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */
457   ierr = PetscOptionsHasName(PETSC_NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr);
458   PetscOptionsEnd();
459   *F = B;
460   PetscFunctionReturn(0);
461 }
462 EXTERN_C_END
463 
464