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: }