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