xref: /petsc/src/mat/impls/aij/seq/umfpack/umfpack.c (revision a093e27361b3b8b9f4c533b65b81b00dcdf12e2c)
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,*ra;
185   UF_long        status;
186   PetscScalar    *av=mat->a;
187 
188   PetscFunctionBegin;
189   if (lu->PetscMatOdering) {
190     ierr = ISGetIndices(r,&ra);CHKERRQ(ierr);
191     ierr = PetscMalloc(m*sizeof(UF_long),&lu->perm_c);CHKERRQ(ierr);
192     /* we cannot simply memcpy on 64 bit archs */
193     for(i = 0; i < m; i++) lu->perm_c[i] = ra[i];
194     ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr);
195   }
196 
197   if(sizeof(UF_long) != sizeof(int)){
198     /* we cannot directly use mat->i and mat->j on 64 bit archs */
199     ierr = PetscMalloc((m+1)*sizeof(UF_long),&lu->ai);CHKERRQ(ierr);
200     ierr = PetscMalloc(mat->nz*sizeof(UF_long),&lu->aj);CHKERRQ(ierr);
201     for(i = 0; i < m + 1; i++) lu->ai[i] = mat->i[i];
202     for(i = 0; i < mat->nz; i++) lu->aj[i] = mat->j[i];
203   }
204   else{
205     lu->ai = (UF_long*)mat->i;
206     lu->aj = (UF_long*)mat->j;
207   }
208 
209   /* print the control parameters */
210 #if defined(PETSC_USE_COMPLEX)
211   if(lu->Control[UMFPACK_PRL] > 1) umfpack_zl_report_control(lu->Control);
212 #else
213   if(lu->Control[UMFPACK_PRL] > 1) umfpack_dl_report_control(lu->Control);
214 #endif
215 
216   /* symbolic factorization of A' */
217   /* ---------------------------------------------------------------------- */
218 #if defined(PETSC_USE_COMPLEX)
219   if (lu->PetscMatOdering) { /* use Petsc row ordering */
220     status = umfpack_zl_qsymbolic(n,m,lu->ai,lu->aj,(double*)av,NULL,
221 				  lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
222   } else { /* use Umfpack col ordering */
223     status = umfpack_zl_symbolic(n,m,lu->ai,lu->aj,(double*)av,NULL,
224 				 &lu->Symbolic,lu->Control,lu->Info);
225   }
226   if (status < 0){
227     umfpack_zl_report_info(lu->Control, lu->Info);
228     umfpack_zl_report_status(lu->Control, status);
229     SETERRQ(PETSC_ERR_LIB,"umfpack_dl_symbolic failed");
230   }
231   /* report sumbolic factorization of A' when Control[PRL] > 3 */
232   (void) umfpack_zl_report_symbolic(lu->Symbolic, lu->Control);
233 #else
234   if (lu->PetscMatOdering) { /* use Petsc row ordering */
235     status = umfpack_dl_qsymbolic(n,m,lu->ai,lu->aj,av,
236 				  lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
237   } else { /* use Umfpack col ordering */
238     status = umfpack_dl_symbolic(n,m,lu->ai,lu->aj,av,
239 				 &lu->Symbolic,lu->Control,lu->Info);
240   }
241   if (status < 0){
242     umfpack_dl_report_info(lu->Control, lu->Info);
243     umfpack_dl_report_status(lu->Control, status);
244     SETERRQ(PETSC_ERR_LIB,"umfpack_dl_symbolic failed");
245   }
246   /* report sumbolic factorization of A' when Control[PRL] > 3 */
247   (void) umfpack_dl_report_symbolic(lu->Symbolic, lu->Control);
248 #endif
249 
250   lu->flg = DIFFERENT_NONZERO_PATTERN;
251   lu->av  = av;
252   lu->CleanUpUMFPACK = PETSC_TRUE;
253   (F)->ops->lufactornumeric  = MatLUFactorNumeric_UMFPACK;
254   PetscFunctionReturn(0);
255 }
256 
257 #undef __FUNCT__
258 #define __FUNCT__ "MatFactorInfo_UMFPACK"
259 PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer)
260 {
261   Mat_UMFPACK    *lu= (Mat_UMFPACK*)A->spptr;
262   PetscErrorCode ierr;
263 
264   PetscFunctionBegin;
265   /* check if matrix is UMFPACK type */
266   if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0);
267 
268   ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr);
269   /* Control parameters used by reporting routiones */
270   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr);
271 
272   /* Control parameters used by symbolic factorization */
273   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);CHKERRQ(ierr);
274   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr);
275   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr);
276   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);CHKERRQ(ierr);
277   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr);
278   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_2BY2_TOLERANCE]: %g\n",lu->Control[UMFPACK_2BY2_TOLERANCE]);CHKERRQ(ierr);
279   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);CHKERRQ(ierr);
280   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);CHKERRQ(ierr);
281 
282   /* Control parameters used by numeric factorization */
283   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr);
284   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);CHKERRQ(ierr);
285   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);CHKERRQ(ierr);
286   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr);
287   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);CHKERRQ(ierr);
288 
289   /* Control parameters used by solve */
290   ierr = PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr);
291 
292   /* mat ordering */
293   if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer,"  UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr);
294 
295   PetscFunctionReturn(0);
296 }
297 
298 #undef __FUNCT__
299 #define __FUNCT__ "MatView_UMFPACK"
300 PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer)
301 {
302   PetscErrorCode    ierr;
303   PetscTruth        iascii;
304   PetscViewerFormat format;
305 
306   PetscFunctionBegin;
307   ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr);
308 
309   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
310   if (iascii) {
311     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
312     if (format == PETSC_VIEWER_ASCII_INFO) {
313       ierr = MatFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr);
314     }
315   }
316   PetscFunctionReturn(0);
317 }
318 
319 /*MC
320   MAT_SOLVER_UMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices
321   via the external package UMFPACK.
322 
323   config/configure.py --download-umfpack to install PETSc to use UMFPACK
324 
325   Consult UMFPACK documentation for more information about the Control parameters
326   which correspond to the options database keys below.
327 
328   Options Database Keys:
329 + -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL]
330 . -mat_umfpack_strategy <AUTO> (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2
331 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
332 . -mat_umfpack_dense_row <0.2>: Control[UMFPACK_DENSE_ROW]
333 . -mat_umfpack_amd_dense <10>: Control[UMFPACK_AMD_DENSE]
334 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
335 . -mat_umfpack_2by2_tolerance <0.01>: Control[UMFPACK_2BY2_TOLERANCE]
336 . -mat_umfpack_fixq <0>: Control[UMFPACK_FIXQ]
337 . -mat_umfpack_aggressive <1>: Control[UMFPACK_AGGRESSIVE]
338 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
339 .  -mat_umfpack_sym_pivot_tolerance <0.001>: Control[UMFPACK_SYM_PIVOT_TOLERANCE]
340 .  -mat_umfpack_scale <NONE> (choose one of) NONE SUM MAX
341 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
342 .  -mat_umfpack_droptol <0>: Control[UMFPACK_DROPTOL]
343 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]
344 
345    Level: beginner
346 
347 .seealso: PCLU, MAT_SOLVER_SUPERLU, MAT_SOLVER_MUMPS, MAT_SOLVER_SPOOLES
348 M*/
349 EXTERN_C_BEGIN
350 #undef __FUNCT__
351 #define __FUNCT__ "MatGetFactor_seqaij_umfpack"
352 PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F)
353 {
354   Mat            B;
355   Mat_UMFPACK    *lu;
356   PetscErrorCode ierr;
357   int            m=A->rmap->n,n=A->cmap->n,idx;
358 
359   const char     *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC","2BY2"},
360                  *scale[]={"NONE","SUM","MAX"};
361   PetscTruth     flg;
362 
363   PetscFunctionBegin;
364   /* Create the factorization matrix F */
365   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
366   ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
367   ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
368   ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
369   ierr = PetscNewLog(B,Mat_UMFPACK,&lu);CHKERRQ(ierr);
370   B->spptr                 = lu;
371   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
372   B->ops->destroy          = MatDestroy_UMFPACK;
373   B->ops->view             = MatView_UMFPACK;
374   B->factor                = MAT_FACTOR_LU;
375   B->assembled             = PETSC_TRUE;  /* required by -ksp_view */
376   B->preallocated          = PETSC_TRUE;
377 
378   /* initializations */
379   /* ------------------------------------------------*/
380   /* get the default control parameters */
381 #if defined(PETSC_USE_COMPLEX)
382   umfpack_zl_defaults(lu->Control);
383 #else
384   umfpack_dl_defaults(lu->Control);
385 #endif
386   lu->perm_c = PETSC_NULL;  /* use defaul UMFPACK col permutation */
387   lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */
388 
389   ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr);
390   /* Control parameters used by reporting routiones */
391   ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr);
392 
393   /* Control parameters for symbolic factorization */
394   ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,4,strategy[0],&idx,&flg);CHKERRQ(ierr);
395   if (flg) {
396     switch (idx){
397     case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break;
398     case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break;
399     case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break;
400     case 3: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2; break;
401     }
402   }
403   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);
404   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);
405   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);
406   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);
407   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);
408   ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);CHKERRQ(ierr);
409   ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);CHKERRQ(ierr);
410 
411   /* Control parameters used by numeric factorization */
412   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);
413   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);
414   ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr);
415   if (flg) {
416     switch (idx){
417     case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
418     case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
419     case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break;
420     }
421   }
422   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);
423   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);
424   ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],PETSC_NULL);CHKERRQ(ierr);
425 
426   /* Control parameters used by solve */
427   ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr);
428 
429   /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */
430   ierr = PetscOptionsHasName(PETSC_NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr);
431   PetscOptionsEnd();
432   *F = B;
433   PetscFunctionReturn(0);
434 }
435 EXTERN_C_END
436 
437