Actual source code: umfpack.c

petsc-3.3-p7 2013-05-11
  2: /* 
  3:    Provides an interface to the UMFPACKv5.1 sparse solver

  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.

 10: */
 11: #include <../src/mat/impls/aij/seq/aij.h>

 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

 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

 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

 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


 75: #define UF_long long long
 76: #define UF_long_max LONG_LONG_MAX
 77: #define UF_long_id "%lld"

 79: EXTERN_C_BEGIN
 80: #include <umfpack.h>
 81: EXTERN_C_END

 83: static const char *const UmfpackOrderingTypes[] = {"CHOLMOD","AMD","GIVEN","METIS","BEST","NONE","USER","UmfpackOrderingTypes","UMFPACK_ORDERING_",0};

 85: typedef struct {
 86:   void         *Symbolic, *Numeric;
 87:   double       Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W;
 88:   PetscInt     *Wi,*perm_c;
 89:   Mat          A;               /* Matrix used for factorization */
 90:   MatStructure flg;
 91:   PetscBool    PetscMatOrdering;

 93:   /* Flag to clean up UMFPACK objects during Destroy */
 94:   PetscBool  CleanUpUMFPACK;
 95: } Mat_UMFPACK;

 99: static PetscErrorCode MatDestroy_UMFPACK(Mat A)
100: {
102:   Mat_UMFPACK    *lu=(Mat_UMFPACK*)A->spptr;

105:   if (lu && lu->CleanUpUMFPACK) {
106:     umfpack_UMF_free_symbolic(&lu->Symbolic);
107:     umfpack_UMF_free_numeric(&lu->Numeric);
108:     PetscFree(lu->Wi);
109:     PetscFree(lu->W);
110:     PetscFree(lu->perm_c);
111:   }
112:   MatDestroy(&lu->A);
113:   PetscFree(A->spptr);
114:   MatDestroy_SeqAIJ(A);
115:   return(0);
116: }

120: static PetscErrorCode MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag)
121: {
122:   Mat_UMFPACK    *lu = (Mat_UMFPACK*)A->spptr;
123:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)lu->A->data;
124:   PetscScalar    *av = a->a,*ba,*xa;
126:   PetscInt       *ai = a->i,*aj = a->j,status;

129:   /* solve Ax = b by umfpack_*_wsolve */
130:   /* ----------------------------------*/

132:   if (!lu->Wi) {  /* first time, allocate working space for wsolve */
133:     PetscMalloc(A->rmap->n*sizeof(PetscInt),&lu->Wi);
134:     PetscMalloc(5*A->rmap->n*sizeof(PetscScalar),&lu->W);
135:   }

137:   VecGetArray(b,&ba);
138:   VecGetArray(x,&xa);
139: #if defined(PETSC_USE_COMPLEX)
140:   status = umfpack_UMF_wsolve(uflag,ai,aj,(PetscReal*)av,NULL,(PetscReal*)xa,NULL,(PetscReal*)ba,NULL,
141:                               lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
142: #else  
143:   status = umfpack_UMF_wsolve(uflag,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
144: #endif
145:   umfpack_UMF_report_info(lu->Control, lu->Info);
146:   if (status < 0){
147:     umfpack_UMF_report_status(lu->Control, status);
148:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_wsolve failed");
149:   }

151:   VecRestoreArray(b,&ba);
152:   VecRestoreArray(x,&xa);
153:   return(0);
154: }

158: static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x)
159: {

163:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
164:   MatSolve_UMFPACK_Private(A,b,x,UMFPACK_Aat);
165:   return(0);
166: }

170: static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x)
171: {

175:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
176:   MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A);
177:   return(0);
178: }

182: static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info)
183: {
184:   Mat_UMFPACK    *lu = (Mat_UMFPACK*)(F)->spptr;
185:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ*)A->data;
186:   PetscInt       *ai = a->i,*aj=a->j,status;
187:   PetscScalar    *av = a->a;

191:   /* numeric factorization of A' */
192:   /* ----------------------------*/

194:   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric){
195:     umfpack_UMF_free_numeric(&lu->Numeric);
196:   }
197: #if defined(PETSC_USE_COMPLEX)
198:   status = umfpack_UMF_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
199: #else
200:   status = umfpack_UMF_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
201: #endif
202:   if (status < 0) {
203:     umfpack_UMF_report_status(lu->Control, status);
204:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_numeric failed");
205:   }
206:   /* report numeric factorization of A' when Control[PRL] > 3 */
207:   (void) umfpack_UMF_report_numeric(lu->Numeric, lu->Control);

209:   PetscObjectReference((PetscObject)A);
210:   MatDestroy(&lu->A);
211:   lu->A = A;
212:   lu->flg = SAME_NONZERO_PATTERN;
213:   lu->CleanUpUMFPACK = PETSC_TRUE;
214:   F->ops->solve          = MatSolve_UMFPACK;
215:   F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
216:   return(0);
217: }

219: /*
220:    Note the r permutation is ignored
221: */
224: static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
225: {
226:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
227:   Mat_UMFPACK    *lu = (Mat_UMFPACK*)(F->spptr);
229:   PetscInt       i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n;
230:   PetscScalar    *av = a->a;
231:   const PetscInt *ra;
232:   PetscInt       status;

235:   if (lu->PetscMatOrdering) {
236:     ISGetIndices(r,&ra);
237:     PetscMalloc(m*sizeof(PetscInt),&lu->perm_c);
238:     /* we cannot simply memcpy on 64 bit archs */
239:     for(i = 0; i < m; i++) lu->perm_c[i] = ra[i];
240:     ISRestoreIndices(r,&ra);
241:   }

243:   /* print the control parameters */
244:   if(lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control);

246:   /* symbolic factorization of A' */
247:   /* ---------------------------------------------------------------------- */
248:   if (lu->PetscMatOrdering) { /* use Petsc row ordering */
249: #if !defined(PETSC_USE_COMPLEX)
250:     status = umfpack_UMF_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
251: #else
252:     status = umfpack_UMF_qsymbolic(n,m,ai,aj,NULL,NULL,
253:                                    lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
254: #endif
255:   } else { /* use Umfpack col ordering */
256: #if !defined(PETSC_USE_COMPLEX)
257:     status = umfpack_UMF_symbolic(n,m,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info);
258: #else
259:     status = umfpack_UMF_symbolic(n,m,ai,aj,NULL,NULL,&lu->Symbolic,lu->Control,lu->Info);
260: #endif
261:   }
262:   if (status < 0){
263:     umfpack_UMF_report_info(lu->Control, lu->Info);
264:     umfpack_UMF_report_status(lu->Control, status);
265:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_symbolic failed");
266:   }
267:   /* report sumbolic factorization of A' when Control[PRL] > 3 */
268:   (void) umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control);

270:   lu->flg = DIFFERENT_NONZERO_PATTERN;
271:   lu->CleanUpUMFPACK = PETSC_TRUE;
272:   (F)->ops->lufactornumeric  = MatLUFactorNumeric_UMFPACK;
273:   return(0);
274: }

278: static PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer)
279: {
280:   Mat_UMFPACK    *lu= (Mat_UMFPACK*)A->spptr;

284:   /* check if matrix is UMFPACK type */
285:   if (A->ops->solve != MatSolve_UMFPACK) return(0);

287:   PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");
288:   /* Control parameters used by reporting routiones */
289:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);

291:   /* Control parameters used by symbolic factorization */
292:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);
293:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);
294:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);
295:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);
296:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);
297:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);
298:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);

300:   /* Control parameters used by numeric factorization */
301:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);
302:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);
303:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);
304:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);
305:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);

307:   /* Control parameters used by solve */
308:   PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);

310:   /* mat ordering */
311:   if (!lu->PetscMatOrdering) {
312:     PetscViewerASCIIPrintf(viewer,"  Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n",UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]]);
313:   }

315:   return(0);
316: }

320: static PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer)
321: {
322:   PetscErrorCode    ierr;
323:   PetscBool         iascii;
324:   PetscViewerFormat format;

327:   MatView_SeqAIJ(A,viewer);

329:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
330:   if (iascii) {
331:     PetscViewerGetFormat(viewer,&format);
332:     if (format == PETSC_VIEWER_ASCII_INFO) {
333:       MatFactorInfo_UMFPACK(A,viewer);
334:     }
335:   }
336:   return(0);
337: }

339: EXTERN_C_BEGIN
342: PetscErrorCode MatFactorGetSolverPackage_seqaij_umfpack(Mat A,const MatSolverPackage *type)
343: {
345:   *type = MATSOLVERUMFPACK;
346:   return(0);
347: }
348: EXTERN_C_END
349: 

351: /*MC
352:   MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices
353:   via the external package UMFPACK.

355:   ./configure --download-umfpack to install PETSc to use UMFPACK

357:   Consult UMFPACK documentation for more information about the Control parameters
358:   which correspond to the options database keys below.

360:   Options Database Keys:
361: + -mat_umfpack_prl                     - UMFPACK print level: Control[UMFPACK_PRL]
362: . -mat_umfpack_strategy <AUTO>         - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2
363: . -mat_umfpack_dense_col <alpha_c>     - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
364: . -mat_umfpack_dense_row <0.2>         - Control[UMFPACK_DENSE_ROW]
365: . -mat_umfpack_amd_dense <10>          - Control[UMFPACK_AMD_DENSE]
366: . -mat_umfpack_block_size <bs>         - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
367: . -mat_umfpack_2by2_tolerance <0.01>   - Control[UMFPACK_2BY2_TOLERANCE]
368: . -mat_umfpack_fixq <0>                - Control[UMFPACK_FIXQ]
369: . -mat_umfpack_aggressive <1>          - Control[UMFPACK_AGGRESSIVE]
370: . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
371: .  -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE]
372: .  -mat_umfpack_scale <NONE>           - (choose one of) NONE SUM MAX
373: . -mat_umfpack_alloc_init <delta>      - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
374: .  -mat_umfpack_droptol <0>            - Control[UMFPACK_DROPTOL]
375: - -mat_umfpack_irstep <maxit>          - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]

377:    Level: beginner

379: .seealso: PCLU, MATSOLVERSUPERLU, MATSOLVERMUMPS, MATSOLVERSPOOLES, PCFactorSetMatSolverPackage(), MatSolverPackage
380: M*/
381: EXTERN_C_BEGIN
384: PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F)
385: {
386:   Mat            B;
387:   Mat_UMFPACK    *lu;
389:   PetscInt       m=A->rmap->n,n=A->cmap->n,idx;

391:   const char     *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC"},
392:                  *scale[]={"NONE","SUM","MAX"};
393:   PetscBool      flg;
394: 
396:   /* Create the factorization matrix F */
397:   MatCreate(((PetscObject)A)->comm,&B);
398:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
399:   MatSetType(B,((PetscObject)A)->type_name);
400:   MatSeqAIJSetPreallocation(B,0,PETSC_NULL);
401:   PetscNewLog(B,Mat_UMFPACK,&lu);
402:   B->spptr                 = lu;
403:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
404:   B->ops->destroy          = MatDestroy_UMFPACK;
405:   B->ops->view             = MatView_UMFPACK;
406:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_umfpack",MatFactorGetSolverPackage_seqaij_umfpack);
407:   B->factortype            = MAT_FACTOR_LU;
408:   B->assembled             = PETSC_TRUE;  /* required by -ksp_view */
409:   B->preallocated          = PETSC_TRUE;
410: 
411:   /* initializations */
412:   /* ------------------------------------------------*/
413:   /* get the default control parameters */
414:   umfpack_UMF_defaults(lu->Control);
415:   lu->perm_c = PETSC_NULL;  /* use defaul UMFPACK col permutation */
416:   lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */

418:   PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"UMFPACK Options","Mat");
419:   /* Control parameters used by reporting routiones */
420:   PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);

422:   /* Control parameters for symbolic factorization */
423:   PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,3,strategy[0],&idx,&flg);
424:   if (flg) {
425:     switch (idx){
426:     case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break;
427:     case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break;
428:     case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break;
429:     }
430:   }
431:   PetscOptionsEList("-mat_umfpack_ordering","Internal ordering method","None",UmfpackOrderingTypes,sizeof UmfpackOrderingTypes/sizeof UmfpackOrderingTypes[0],UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]],&idx,&flg);
432:   if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx;
433:   PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],PETSC_NULL);
434:   PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],PETSC_NULL);
435:   PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],PETSC_NULL);
436:   PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],PETSC_NULL);
437:   PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);
438:   PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);

440:   /* Control parameters used by numeric factorization */
441:   PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],PETSC_NULL);
442:   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);
443:   PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);
444:   if (flg) {
445:     switch (idx){
446:     case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
447:     case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
448:     case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break;
449:     }
450:   }
451:   PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);
452:   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);
453:   PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],PETSC_NULL);

455:   /* Control parameters used by solve */
456:   PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);
457: 
458:   /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */
459:   PetscOptionsHasName(PETSC_NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOrdering);
460:   PetscOptionsEnd();
461:   *F = B;
462:   return(0);
463: }
464: EXTERN_C_END