Actual source code: umfpack.c

petsc-master 2015-07-26
Report Typos and Errors
  2: /*
  3:    Provides an interface to the UMFPACK sparse solver available through SuiteSparse version 4.2.1

  5:    When build with PETSC_USE_64BIT_INDICES this will use Suitesparse_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 one cannot use 64BIT_INDICES on 32bit machines [as Suitesparse_long is 32bit only]

 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: /* the type casts are needed because PetscInt is long long while SuiteSparse_long is long and compilers warn even when they are identical */
 18: #define umfpack_UMF_wsolve(a,b,c,d,e,f,g,h,i,j,k,l,m,n) umfpack_zl_wsolve(a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,d,e,f,g,h,i,(SuiteSparse_long*)j,k,l,(SuiteSparse_long*)m,n)
 19: #define umfpack_UMF_numeric(a,b,c,d,e,f,g,h)          umfpack_zl_numeric((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e,f,g,h)
 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(a,b,c,d,e,f,g,h,i,j)    umfpack_zl_qsymbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,(SuiteSparse_long*)g,h,i,j)
 26: #define umfpack_UMF_symbolic(a,b,c,d,e,f,g,h,i)       umfpack_zl_symbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,g,h,i)
 27: #define umfpack_UMF_defaults                          umfpack_zl_defaults

 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(a,b,c,d,e,f,g,h,i,j,k)  umfpack_dl_wsolve(a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,d,e,f,g,h,i,(SuiteSparse_long*)j,k)
 33: #define umfpack_UMF_numeric(a,b,c,d,e,f,g)         umfpack_dl_numeric((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e,f,g)
 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(a,b,c,d,e,f,g,h,i)   umfpack_dl_qsymbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,(SuiteSparse_long*)f,g,h,i)
 40: #define umfpack_UMF_symbolic(a,b,c,d,e,f,g,h)      umfpack_dl_symbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,g,h)
 41: #define umfpack_UMF_defaults                       umfpack_dl_defaults
 42: #endif

 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

 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

 75: EXTERN_C_BEGIN
 76: #include <umfpack.h>
 77: EXTERN_C_END

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

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

 89:   /* Flag to clean up UMFPACK objects during Destroy */
 90:   PetscBool CleanUpUMFPACK;
 91: } Mat_UMFPACK;

 95: static PetscErrorCode MatDestroy_UMFPACK(Mat A)
 96: {
 98:   Mat_UMFPACK    *lu=(Mat_UMFPACK*)A->spptr;

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

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

126:   /* solve Ax = b by umfpack_*_wsolve */
127:   /* ----------------------------------*/

129:   if (!lu->Wi) {  /* first time, allocate working space for wsolve */
130:     PetscMalloc1(A->rmap->n,&lu->Wi);
131:     PetscMalloc1(5*A->rmap->n,&lu->W);
132:   }

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

147:   VecRestoreArrayRead(b,&ba);
148:   VecRestoreArray(x,&xa);
149:   return(0);
150: }

154: static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x)
155: {

159:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
160:   MatSolve_UMFPACK_Private(A,b,x,UMFPACK_Aat);
161:   return(0);
162: }

166: static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x)
167: {

171:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
172:   MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A);
173:   return(0);
174: }

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

187:   /* numeric factorization of A' */
188:   /* ----------------------------*/

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

205:   PetscObjectReference((PetscObject)A);
206:   MatDestroy(&lu->A);

208:   lu->A                  = A;
209:   lu->flg                = SAME_NONZERO_PATTERN;
210:   lu->CleanUpUMFPACK     = PETSC_TRUE;
211:   F->ops->solve          = MatSolve_UMFPACK;
212:   F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
213:   return(0);
214: }

216: /*
217:    Note the r permutation is ignored
218: */
221: static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
222: {
223:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ*)A->data;
224:   Mat_UMFPACK    *lu = (Mat_UMFPACK*)(F->spptr);
226:   PetscInt       i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n;
227: #if !defined(PETSC_USE_COMPLEX)
228:   PetscScalar    *av = a->a;
229: #endif
230:   const PetscInt *ra;
231:   PetscInt       status;

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

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

245:   /* symbolic factorization of A' */
246:   /* ---------------------------------------------------------------------- */
247:   if (lu->PetscMatOrdering) { /* use Petsc row ordering */
248: #if !defined(PETSC_USE_COMPLEX)
249:     status = umfpack_UMF_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
250: #else
251:     status = umfpack_UMF_qsymbolic(n,m,ai,aj,NULL,NULL,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,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info);
256: #else
257:     status = umfpack_UMF_symbolic(n,m,ai,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_COMM_SELF,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);

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

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

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

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

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

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

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

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

317: static PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer)
318: {
319:   PetscErrorCode    ierr;
320:   PetscBool         iascii;
321:   PetscViewerFormat format;

324:   MatView_SeqAIJ(A,viewer);

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

338: PetscErrorCode MatFactorGetSolverPackage_seqaij_umfpack(Mat A,const MatSolverPackage *type)
339: {
341:   *type = MATSOLVERUMFPACK;
342:   return(0);
343: }


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

350:   Use ./configure --download-suitesparse to install PETSc to use UMFPACK

352:   Use -pc_type lu -pc_factor_mat_solver_package umfpack to us this direct solver

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

357:   Options Database Keys:
358: + -mat_umfpack_ordering                - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE
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]

375:    Level: beginner

377:    Note: UMFPACK is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html

379: .seealso: PCLU, MATSOLVERSUPERLU, MATSOLVERMUMPS, PCFactorSetMatSolverPackage(), MatSolverPackage
380: M*/

384: PETSC_EXTERN 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:   const char *scale[]   ={"NONE","SUM","MAX"};
393:   PetscBool  flg;

396:   /* Create the factorization matrix F */
397:   MatCreate(PetscObjectComm((PetscObject)A),&B);
398:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
399:   MatSetType(B,((PetscObject)A)->type_name);
400:   MatSeqAIJSetPreallocation(B,0,NULL);
401:   PetscNewLog(B,&lu);

403:   B->spptr                 = lu;
404:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
405:   B->ops->destroy          = MatDestroy_UMFPACK;
406:   B->ops->view             = MatView_UMFPACK;
407:   B->ops->matsolve         = NULL;

409:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_umfpack);

411:   B->factortype   = MAT_FACTOR_LU;
412:   B->assembled    = PETSC_TRUE;           /* required by -ksp_view */
413:   B->preallocated = PETSC_TRUE;

415:   /* initializations */
416:   /* ------------------------------------------------*/
417:   /* get the default control parameters */
418:   umfpack_UMF_defaults(lu->Control);
419:   lu->perm_c                  = NULL; /* use defaul UMFPACK col permutation */
420:   lu->Control[UMFPACK_IRSTEP] = 0;          /* max num of iterative refinement steps to attempt */

422:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"UMFPACK Options","Mat");
423:   /* Control parameters used by reporting routiones */
424:   PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],NULL);

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

444:   /* Control parameters used by numeric factorization */
445:   PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],NULL);
446:   PetscOptionsReal("-mat_umfpack_sym_pivot_tolerance","Control[UMFPACK_SYM_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],&lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],NULL);
447:   PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);
448:   if (flg) {
449:     switch (idx) {
450:     case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
451:     case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
452:     case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break;
453:     }
454:   }
455:   PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL);
456:   PetscOptionsReal("-mat_umfpack_front_alloc_init","Control[UMFPACK_FRONT_ALLOC_INIT]","None",lu->Control[UMFPACK_FRONT_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL);
457:   PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],NULL);

459:   /* Control parameters used by solve */
460:   PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],NULL);

462:   /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */
463:   PetscOptionsHasName(NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOrdering);
464:   PetscOptionsEnd();
465:   *F = B;
466:   return(0);
467: }

469: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*);
470: PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*);
471: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat,MatFactorType,Mat*);

475: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_SuiteSparse(void)
476: {

480:   MatSolverPackageRegister(MATSOLVERUMFPACK,MATSEQAIJ,      MAT_FACTOR_LU,MatGetFactor_seqaij_umfpack);
481:   MatSolverPackageRegister(MATSOLVERCHOLMOD,MATSEQAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_seqaij_cholmod);
482:   MatSolverPackageRegister(MATSOLVERCHOLMOD,MATSEQSBAIJ,      MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_cholmod);
483:   MatSolverPackageRegister(MATSOLVERKLU,MATSEQAIJ,          MAT_FACTOR_LU,MatGetFactor_seqaij_klu);
484:   return(0);
485: }