Actual source code: umfpack.c


  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;

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

 92: static PetscErrorCode MatDestroy_UMFPACK(Mat A)
 93: {
 94:   Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;

 96:   if (lu->CleanUpUMFPACK) {
 97:     umfpack_UMF_free_symbolic(&lu->Symbolic);
 98:     umfpack_UMF_free_numeric(&lu->Numeric);
 99:     PetscFree(lu->Wi);
100:     PetscFree(lu->W);
101:     PetscFree(lu->perm_c);
102:   }
103:   MatDestroy(&lu->A);
104:   PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL);
105:   PetscFree(A->data);
106:   return 0;
107: }

109: static PetscErrorCode MatSolve_UMFPACK_Private(Mat A, Vec b, Vec x, int uflag)
110: {
111:   Mat_UMFPACK       *lu = (Mat_UMFPACK *)A->data;
112:   Mat_SeqAIJ        *a  = (Mat_SeqAIJ *)lu->A->data;
113:   PetscScalar       *av = a->a, *xa;
114:   const PetscScalar *ba;
115:   PetscInt          *ai = a->i, *aj = a->j, status;
116:   static PetscBool   cite = PETSC_FALSE;

118:   if (!A->rmap->n) return 0;
119:   PetscCall(PetscCitationsRegister("@article{davis2004algorithm,\n  title={Algorithm 832: {UMFPACK} V4.3---An Unsymmetric-Pattern Multifrontal Method},\n  author={Davis, Timothy A},\n  journal={ACM Transactions on Mathematical Software (TOMS)},\n  "
120:                                    "volume={30},\n  number={2},\n  pages={196--199},\n  year={2004},\n  publisher={ACM}\n}\n",
121:                                    &cite));
122:   /* solve Ax = b by umfpack_*_wsolve */
123:   /* ----------------------------------*/

125:   if (!lu->Wi) { /* first time, allocate working space for wsolve */
126:     PetscMalloc1(A->rmap->n, &lu->Wi);
127:     PetscMalloc1(5 * A->rmap->n, &lu->W);
128:   }

130:   VecGetArrayRead(b, &ba);
131:   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, lu->Numeric, lu->Control, lu->Info, lu->Wi, lu->W);
134: #else
135:   status = umfpack_UMF_wsolve(uflag, ai, aj, av, xa, ba, lu->Numeric, lu->Control, lu->Info, lu->Wi, lu->W);
136: #endif
137:   umfpack_UMF_report_info(lu->Control, lu->Info);
138:   if (status < 0) {
139:     umfpack_UMF_report_status(lu->Control, status);
140:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_wsolve failed");
141:   }

143:   VecRestoreArrayRead(b, &ba);
144:   VecRestoreArray(x, &xa);
145:   return 0;
146: }

148: static PetscErrorCode MatSolve_UMFPACK(Mat A, Vec b, Vec x)
149: {
150:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
151:   MatSolve_UMFPACK_Private(A, b, x, UMFPACK_Aat);
152:   return 0;
153: }

155: static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A, Vec b, Vec x)
156: {
157:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
158:   MatSolve_UMFPACK_Private(A, b, x, UMFPACK_A);
159:   return 0;
160: }

162: static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F, Mat A, const MatFactorInfo *info)
163: {
164:   Mat_UMFPACK *lu = (Mat_UMFPACK *)(F)->data;
165:   Mat_SeqAIJ  *a  = (Mat_SeqAIJ *)A->data;
166:   PetscInt    *ai = a->i, *aj = a->j, status;
167:   PetscScalar *av = a->a;

169:   if (!A->rmap->n) return 0;
170:   /* numeric factorization of A' */
171:   /* ----------------------------*/

173:   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) umfpack_UMF_free_numeric(&lu->Numeric);
174: #if defined(PETSC_USE_COMPLEX)
175:   status = umfpack_UMF_numeric(ai, aj, (double *)av, NULL, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info);
176: #else
177:   status = umfpack_UMF_numeric(ai, aj, av, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info);
178: #endif
179:   if (status < 0) {
180:     umfpack_UMF_report_status(lu->Control, status);
181:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_numeric failed");
182:   }
183:   /* report numeric factorization of A' when Control[PRL] > 3 */
184:   (void)umfpack_UMF_report_numeric(lu->Numeric, lu->Control);

186:   PetscObjectReference((PetscObject)A);
187:   MatDestroy(&lu->A);

189:   lu->A                  = A;
190:   lu->flg                = SAME_NONZERO_PATTERN;
191:   lu->CleanUpUMFPACK     = PETSC_TRUE;
192:   F->ops->solve          = MatSolve_UMFPACK;
193:   F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
194:   return 0;
195: }

197: static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
198: {
199:   Mat_SeqAIJ  *a  = (Mat_SeqAIJ *)A->data;
200:   Mat_UMFPACK *lu = (Mat_UMFPACK *)(F->data);
201:   PetscInt     i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n, status, idx;
202: #if !defined(PETSC_USE_COMPLEX)
203:   PetscScalar *av = a->a;
204: #endif
205:   const PetscInt *ra;
206:   const char     *strategy[] = {"AUTO", "UNSYMMETRIC", "SYMMETRIC"};
207:   const char     *scale[]    = {"NONE", "SUM", "MAX"};
208:   PetscBool       flg;

210:   (F)->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK;
211:   if (!n) return 0;

213:   /* Set options to F */
214:   PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "UMFPACK Options", "Mat");
215:   /* Control parameters used by reporting routiones */
216:   PetscOptionsReal("-mat_umfpack_prl", "Control[UMFPACK_PRL]", "None", lu->Control[UMFPACK_PRL], &lu->Control[UMFPACK_PRL], NULL);

218:   /* Control parameters for symbolic factorization */
219:   PetscOptionsEList("-mat_umfpack_strategy", "ordering and pivoting strategy", "None", strategy, 3, strategy[0], &idx, &flg);
220:   if (flg) {
221:     switch (idx) {
222:     case 0:
223:       lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO;
224:       break;
225:     case 1:
226:       lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC;
227:       break;
228:     case 2:
229:       lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC;
230:       break;
231:     }
232:   }
233:   PetscOptionsEList("-mat_umfpack_ordering", "Internal ordering method", "None", UmfpackOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(UmfpackOrderingTypes), UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]], &idx, &flg);
234:   if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx;
235:   PetscOptionsReal("-mat_umfpack_dense_col", "Control[UMFPACK_DENSE_COL]", "None", lu->Control[UMFPACK_DENSE_COL], &lu->Control[UMFPACK_DENSE_COL], NULL);
236:   PetscOptionsReal("-mat_umfpack_dense_row", "Control[UMFPACK_DENSE_ROW]", "None", lu->Control[UMFPACK_DENSE_ROW], &lu->Control[UMFPACK_DENSE_ROW], NULL);
237:   PetscOptionsReal("-mat_umfpack_amd_dense", "Control[UMFPACK_AMD_DENSE]", "None", lu->Control[UMFPACK_AMD_DENSE], &lu->Control[UMFPACK_AMD_DENSE], NULL);
238:   PetscOptionsReal("-mat_umfpack_block_size", "Control[UMFPACK_BLOCK_SIZE]", "None", lu->Control[UMFPACK_BLOCK_SIZE], &lu->Control[UMFPACK_BLOCK_SIZE], NULL);
239:   PetscOptionsReal("-mat_umfpack_fixq", "Control[UMFPACK_FIXQ]", "None", lu->Control[UMFPACK_FIXQ], &lu->Control[UMFPACK_FIXQ], NULL);
240:   PetscOptionsReal("-mat_umfpack_aggressive", "Control[UMFPACK_AGGRESSIVE]", "None", lu->Control[UMFPACK_AGGRESSIVE], &lu->Control[UMFPACK_AGGRESSIVE], NULL);

242:   /* Control parameters used by numeric factorization */
243:   PetscOptionsReal("-mat_umfpack_pivot_tolerance", "Control[UMFPACK_PIVOT_TOLERANCE]", "None", lu->Control[UMFPACK_PIVOT_TOLERANCE], &lu->Control[UMFPACK_PIVOT_TOLERANCE], NULL);
244:   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);
245:   PetscOptionsEList("-mat_umfpack_scale", "Control[UMFPACK_SCALE]", "None", scale, 3, scale[0], &idx, &flg);
246:   if (flg) {
247:     switch (idx) {
248:     case 0:
249:       lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE;
250:       break;
251:     case 1:
252:       lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM;
253:       break;
254:     case 2:
255:       lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX;
256:       break;
257:     }
258:   }
259:   PetscOptionsReal("-mat_umfpack_alloc_init", "Control[UMFPACK_ALLOC_INIT]", "None", lu->Control[UMFPACK_ALLOC_INIT], &lu->Control[UMFPACK_ALLOC_INIT], NULL);
260:   PetscOptionsReal("-mat_umfpack_front_alloc_init", "Control[UMFPACK_FRONT_ALLOC_INIT]", "None", lu->Control[UMFPACK_FRONT_ALLOC_INIT], &lu->Control[UMFPACK_ALLOC_INIT], NULL);
261:   PetscOptionsReal("-mat_umfpack_droptol", "Control[UMFPACK_DROPTOL]", "None", lu->Control[UMFPACK_DROPTOL], &lu->Control[UMFPACK_DROPTOL], NULL);

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

267:   if (r) {
268:     ISGetIndices(r, &ra);
269:     PetscMalloc1(m, &lu->perm_c);
270:     /* we cannot simply memcpy on 64 bit archs */
271:     for (i = 0; i < m; i++) lu->perm_c[i] = ra[i];
272:     ISRestoreIndices(r, &ra);
273:   }

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

278:   /* symbolic factorization of A' */
279:   /* ---------------------------------------------------------------------- */
280:   if (r) { /* use Petsc row ordering */
281: #if !defined(PETSC_USE_COMPLEX)
282:     status = umfpack_UMF_qsymbolic(n, m, ai, aj, av, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info);
283: #else
284:     status = umfpack_UMF_qsymbolic(n, m, ai, aj, NULL, NULL, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info);
285: #endif
286:   } else { /* use Umfpack col ordering */
287: #if !defined(PETSC_USE_COMPLEX)
288:     status = umfpack_UMF_symbolic(n, m, ai, aj, av, &lu->Symbolic, lu->Control, lu->Info);
289: #else
290:     status = umfpack_UMF_symbolic(n, m, ai, aj, NULL, NULL, &lu->Symbolic, lu->Control, lu->Info);
291: #endif
292:   }
293:   if (status < 0) {
294:     umfpack_UMF_report_info(lu->Control, lu->Info);
295:     umfpack_UMF_report_status(lu->Control, status);
296:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_symbolic failed");
297:   }
298:   /* report sumbolic factorization of A' when Control[PRL] > 3 */
299:   (void)umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control);

301:   lu->flg            = DIFFERENT_NONZERO_PATTERN;
302:   lu->CleanUpUMFPACK = PETSC_TRUE;
303:   return 0;
304: }

306: static PetscErrorCode MatView_Info_UMFPACK(Mat A, PetscViewer viewer)
307: {
308:   Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;

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

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

317:   /* Control parameters used by symbolic factorization */
318:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_STRATEGY]: %g\n", lu->Control[UMFPACK_STRATEGY]);
319:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_DENSE_COL]: %g\n", lu->Control[UMFPACK_DENSE_COL]);
320:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_DENSE_ROW]: %g\n", lu->Control[UMFPACK_DENSE_ROW]);
321:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_AMD_DENSE]: %g\n", lu->Control[UMFPACK_AMD_DENSE]);
322:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_BLOCK_SIZE]: %g\n", lu->Control[UMFPACK_BLOCK_SIZE]);
323:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_FIXQ]: %g\n", lu->Control[UMFPACK_FIXQ]);
324:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_AGGRESSIVE]: %g\n", lu->Control[UMFPACK_AGGRESSIVE]);

326:   /* Control parameters used by numeric factorization */
327:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_PIVOT_TOLERANCE]);
328:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);
329:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_SCALE]: %g\n", lu->Control[UMFPACK_SCALE]);
330:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_ALLOC_INIT]: %g\n", lu->Control[UMFPACK_ALLOC_INIT]);
331:   PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_DROPTOL]: %g\n", lu->Control[UMFPACK_DROPTOL]);

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

336:   /* mat ordering */
337:   if (!lu->perm_c) PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n", UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]]);
338:   return 0;
339: }

341: static PetscErrorCode MatView_UMFPACK(Mat A, PetscViewer viewer)
342: {
343:   PetscBool         iascii;
344:   PetscViewerFormat format;

346:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
347:   if (iascii) {
348:     PetscViewerGetFormat(viewer, &format);
349:     if (format == PETSC_VIEWER_ASCII_INFO) MatView_Info_UMFPACK(A, viewer);
350:   }
351:   return 0;
352: }

354: PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A, MatSolverType *type)
355: {
356:   *type = MATSOLVERUMFPACK;
357:   return 0;
358: }

360: /*MC
361:   MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers, LU, for sequential matrices
362:   via the external package UMFPACK.

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

366:   Use -pc_type lu -pc_factor_mat_solver_type umfpack to use this direct solver

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

371:   Options Database Keys:
372: + -mat_umfpack_ordering                - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE
373: . -mat_umfpack_prl                     - UMFPACK print level: Control[UMFPACK_PRL]
374: . -mat_umfpack_strategy <AUTO>         - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2
375: . -mat_umfpack_dense_col <alpha_c>     - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
376: . -mat_umfpack_dense_row <0.2>         - Control[UMFPACK_DENSE_ROW]
377: . -mat_umfpack_amd_dense <10>          - Control[UMFPACK_AMD_DENSE]
378: . -mat_umfpack_block_size <bs>         - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
379: . -mat_umfpack_2by2_tolerance <0.01>   - Control[UMFPACK_2BY2_TOLERANCE]
380: . -mat_umfpack_fixq <0>                - Control[UMFPACK_FIXQ]
381: . -mat_umfpack_aggressive <1>          - Control[UMFPACK_AGGRESSIVE]
382: . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
383: . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE]
384: . -mat_umfpack_scale <NONE>           - (choose one of) NONE SUM MAX
385: . -mat_umfpack_alloc_init <delta>      - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
386: . -mat_umfpack_droptol <0>            - Control[UMFPACK_DROPTOL]
387: - -mat_umfpack_irstep <maxit>          - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]

389:    Level: beginner

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

393: .seealso: `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`
394: M*/

396: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A, MatFactorType ftype, Mat *F)
397: {
398:   Mat          B;
399:   Mat_UMFPACK *lu;
400:   PetscInt     m = A->rmap->n, n = A->cmap->n;

402:   /* Create the factorization matrix F */
403:   MatCreate(PetscObjectComm((PetscObject)A), &B);
404:   MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n);
405:   PetscStrallocpy("umfpack", &((PetscObject)B)->type_name);
406:   MatSetUp(B);

408:   PetscNew(&lu);

410:   B->data                  = lu;
411:   B->ops->getinfo          = MatGetInfo_External;
412:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
413:   B->ops->destroy          = MatDestroy_UMFPACK;
414:   B->ops->view             = MatView_UMFPACK;
415:   B->ops->matsolve         = NULL;

417:   PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_umfpack);

419:   B->factortype   = MAT_FACTOR_LU;
420:   B->assembled    = PETSC_TRUE; /* required by -ksp_view */
421:   B->preallocated = PETSC_TRUE;

423:   PetscFree(B->solvertype);
424:   PetscStrallocpy(MATSOLVERUMFPACK, &B->solvertype);
425:   B->canuseordering = PETSC_TRUE;
426:   PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]);

428:   /* initializations */
429:   /* ------------------------------------------------*/
430:   /* get the default control parameters */
431:   umfpack_UMF_defaults(lu->Control);
432:   lu->perm_c                  = NULL; /* use defaul UMFPACK col permutation */
433:   lu->Control[UMFPACK_IRSTEP] = 0;    /* max num of iterative refinement steps to attempt */

435:   *F = B;
436:   return 0;
437: }

439: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat, MatFactorType, Mat *);
440: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat, MatFactorType, Mat *);
441: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat, MatFactorType, Mat *);
442: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat, MatFactorType, Mat *);

444: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void)
445: {
446:   MatSolverTypeRegister(MATSOLVERUMFPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_umfpack);
447:   MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqaij_cholmod);
448:   MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_cholmod);
449:   MatSolverTypeRegister(MATSOLVERKLU, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_klu);
450:   MatSolverTypeRegister(MATSOLVERSPQR, MATSEQAIJ, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr);
451:   if (!PetscDefined(USE_COMPLEX)) MatSolverTypeRegister(MATSOLVERSPQR, MATNORMAL, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr);
452:   MatSolverTypeRegister(MATSOLVERSPQR, MATNORMALHERMITIAN, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr);
453:   return 0;
454: }