Actual source code: umfpack.c
1: /*
2: Provides an interface to the UMFPACK sparse solver available through SuiteSparse version 4.2.1
4: When build with PETSC_USE_64BIT_INDICES this will use Suitesparse_long as the
5: integer type in UMFPACK, otherwise it will use int. This means
6: all integers in this file as simply declared as PetscInt. Also it means
7: that one cannot use 64BIT_INDICES on 32-bit pointer systems [as Suitesparse_long is 32-bit only]
9: */
10: #include <../src/mat/impls/aij/seq/aij.h>
12: #if defined(PETSC_USE_64BIT_INDICES)
13: #if defined(PETSC_USE_COMPLEX)
14: #define umfpack_UMF_free_symbolic umfpack_zl_free_symbolic
15: #define umfpack_UMF_free_numeric umfpack_zl_free_numeric
16: /* the type casts are needed because PetscInt is long long while SuiteSparse_long is long and compilers warn even when they are identical */
17: #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)
18: #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)
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(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)
25: #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)
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(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)
32: #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)
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(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)
39: #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)
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
74: EXTERN_C_BEGIN
75: #include <umfpack.h>
76: EXTERN_C_END
78: static const char *const UmfpackOrderingTypes[] = {"CHOLMOD", "AMD", "GIVEN", "METIS", "BEST", "NONE", "USER", "UmfpackOrderingTypes", "UMFPACK_ORDERING_", 0};
80: typedef struct {
81: void *Symbolic, *Numeric;
82: double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL], *W;
83: PetscInt *Wi, *perm_c;
84: Mat A; /* Matrix used for factorization */
85: MatStructure flg;
87: /* Flag to clean up UMFPACK objects during Destroy */
88: PetscBool CleanUpUMFPACK;
89: } Mat_UMFPACK;
91: static PetscErrorCode MatDestroy_UMFPACK(Mat A)
92: {
93: Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;
95: PetscFunctionBegin;
96: if (lu->CleanUpUMFPACK) {
97: umfpack_UMF_free_symbolic(&lu->Symbolic);
98: umfpack_UMF_free_numeric(&lu->Numeric);
99: PetscCall(PetscFree(lu->Wi));
100: PetscCall(PetscFree(lu->W));
101: PetscCall(PetscFree(lu->perm_c));
102: }
103: PetscCall(MatDestroy(&lu->A));
104: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
105: PetscCall(PetscFree(A->data));
106: PetscFunctionReturn(PETSC_SUCCESS);
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;
116: int status;
117: static PetscBool cite = PETSC_FALSE;
119: PetscFunctionBegin;
120: if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS);
121: 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 "
122: "volume={30},\n number={2},\n pages={196--199},\n year={2004},\n publisher={ACM}\n}\n",
123: &cite));
124: /* solve Ax = b by umfpack_*_wsolve */
126: if (!lu->Wi) { /* first time, allocate working space for wsolve */
127: PetscCall(PetscMalloc1(A->rmap->n, &lu->Wi));
128: PetscCall(PetscMalloc1(5 * A->rmap->n, &lu->W));
129: }
131: PetscCall(VecGetArrayRead(b, &ba));
132: PetscCall(VecGetArray(x, &xa));
133: #if defined(PETSC_USE_COMPLEX)
134: 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);
135: #else
136: status = umfpack_UMF_wsolve(uflag, ai, aj, av, xa, ba, lu->Numeric, lu->Control, lu->Info, lu->Wi, lu->W);
137: #endif
138: umfpack_UMF_report_info(lu->Control, lu->Info);
139: if (status < 0) {
140: umfpack_UMF_report_status(lu->Control, status);
141: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_wsolve failed");
142: }
144: PetscCall(VecRestoreArrayRead(b, &ba));
145: PetscCall(VecRestoreArray(x, &xa));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode MatSolve_UMFPACK(Mat A, Vec b, Vec x)
150: {
151: PetscFunctionBegin;
152: /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
153: PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_Aat));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A, Vec b, Vec x)
158: {
159: PetscFunctionBegin;
160: /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
161: PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_A));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F, Mat A, const MatFactorInfo *info)
166: {
167: Mat_UMFPACK *lu = (Mat_UMFPACK *)F->data;
168: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
169: PetscInt *ai = a->i, *aj = a->j;
170: int status;
171: PetscScalar *av = a->a;
173: PetscFunctionBegin;
174: if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS);
175: /* numeric factorization of A' */
177: if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) umfpack_UMF_free_numeric(&lu->Numeric);
178: #if defined(PETSC_USE_COMPLEX)
179: status = umfpack_UMF_numeric(ai, aj, (double *)av, NULL, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info);
180: #else
181: status = umfpack_UMF_numeric(ai, aj, av, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info);
182: #endif
183: if (status < 0) {
184: umfpack_UMF_report_status(lu->Control, status);
185: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_numeric failed");
186: }
187: /* report numeric factorization of A' when Control[PRL] > 3 */
188: (void)umfpack_UMF_report_numeric(lu->Numeric, lu->Control);
190: PetscCall(PetscObjectReference((PetscObject)A));
191: PetscCall(MatDestroy(&lu->A));
193: lu->A = A;
194: lu->flg = SAME_NONZERO_PATTERN;
195: lu->CleanUpUMFPACK = PETSC_TRUE;
196: F->ops->solve = MatSolve_UMFPACK;
197: F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
202: {
203: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
204: Mat_UMFPACK *lu = (Mat_UMFPACK *)F->data;
205: PetscInt i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n, idx;
206: int status;
207: #if !defined(PETSC_USE_COMPLEX)
208: PetscScalar *av = a->a;
209: #endif
210: const PetscInt *ra;
211: const char *strategy[] = {"AUTO", "UNSYMMETRIC", "SYMMETRIC"};
212: const char *scale[] = {"NONE", "SUM", "MAX"};
213: PetscBool flg;
215: PetscFunctionBegin;
216: F->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK;
217: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
219: /* Set options to F */
220: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "UMFPACK Options", "Mat");
221: /* Control parameters used by reporting routiones */
222: PetscCall(PetscOptionsReal("-mat_umfpack_prl", "Control[UMFPACK_PRL]", "None", lu->Control[UMFPACK_PRL], &lu->Control[UMFPACK_PRL], NULL));
224: /* Control parameters for symbolic factorization */
225: PetscCall(PetscOptionsEList("-mat_umfpack_strategy", "ordering and pivoting strategy", "None", strategy, 3, strategy[0], &idx, &flg));
226: if (flg) {
227: switch (idx) {
228: case 0:
229: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO;
230: break;
231: case 1:
232: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC;
233: break;
234: case 2:
235: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC;
236: break;
237: }
238: }
239: PetscCall(PetscOptionsEList("-mat_umfpack_ordering", "Internal ordering method", "None", UmfpackOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(UmfpackOrderingTypes), UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]], &idx, &flg));
240: if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx;
241: PetscCall(PetscOptionsReal("-mat_umfpack_dense_col", "Control[UMFPACK_DENSE_COL]", "None", lu->Control[UMFPACK_DENSE_COL], &lu->Control[UMFPACK_DENSE_COL], NULL));
242: PetscCall(PetscOptionsReal("-mat_umfpack_dense_row", "Control[UMFPACK_DENSE_ROW]", "None", lu->Control[UMFPACK_DENSE_ROW], &lu->Control[UMFPACK_DENSE_ROW], NULL));
243: PetscCall(PetscOptionsReal("-mat_umfpack_amd_dense", "Control[UMFPACK_AMD_DENSE]", "None", lu->Control[UMFPACK_AMD_DENSE], &lu->Control[UMFPACK_AMD_DENSE], NULL));
244: PetscCall(PetscOptionsReal("-mat_umfpack_block_size", "Control[UMFPACK_BLOCK_SIZE]", "None", lu->Control[UMFPACK_BLOCK_SIZE], &lu->Control[UMFPACK_BLOCK_SIZE], NULL));
245: PetscCall(PetscOptionsReal("-mat_umfpack_fixq", "Control[UMFPACK_FIXQ]", "None", lu->Control[UMFPACK_FIXQ], &lu->Control[UMFPACK_FIXQ], NULL));
246: PetscCall(PetscOptionsReal("-mat_umfpack_aggressive", "Control[UMFPACK_AGGRESSIVE]", "None", lu->Control[UMFPACK_AGGRESSIVE], &lu->Control[UMFPACK_AGGRESSIVE], NULL));
248: /* Control parameters used by numeric factorization */
249: PetscCall(PetscOptionsReal("-mat_umfpack_pivot_tolerance", "Control[UMFPACK_PIVOT_TOLERANCE]", "None", lu->Control[UMFPACK_PIVOT_TOLERANCE], &lu->Control[UMFPACK_PIVOT_TOLERANCE], NULL));
250: PetscCall(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));
251: PetscCall(PetscOptionsEList("-mat_umfpack_scale", "Control[UMFPACK_SCALE]", "None", scale, 3, scale[0], &idx, &flg));
252: if (flg) {
253: switch (idx) {
254: case 0:
255: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE;
256: break;
257: case 1:
258: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM;
259: break;
260: case 2:
261: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX;
262: break;
263: }
264: }
265: PetscCall(PetscOptionsReal("-mat_umfpack_alloc_init", "Control[UMFPACK_ALLOC_INIT]", "None", lu->Control[UMFPACK_ALLOC_INIT], &lu->Control[UMFPACK_ALLOC_INIT], NULL));
266: PetscCall(PetscOptionsReal("-mat_umfpack_front_alloc_init", "Control[UMFPACK_FRONT_ALLOC_INIT]", "None", lu->Control[UMFPACK_FRONT_ALLOC_INIT], &lu->Control[UMFPACK_ALLOC_INIT], NULL));
267: PetscCall(PetscOptionsReal("-mat_umfpack_droptol", "Control[UMFPACK_DROPTOL]", "None", lu->Control[UMFPACK_DROPTOL], &lu->Control[UMFPACK_DROPTOL], NULL));
269: /* Control parameters used by solve */
270: PetscCall(PetscOptionsReal("-mat_umfpack_irstep", "Control[UMFPACK_IRSTEP]", "None", lu->Control[UMFPACK_IRSTEP], &lu->Control[UMFPACK_IRSTEP], NULL));
271: PetscOptionsEnd();
273: if (r) {
274: PetscCall(ISGetIndices(r, &ra));
275: PetscCall(PetscMalloc1(m, &lu->perm_c));
276: /* we cannot simply memcpy on 64-bit archs */
277: for (i = 0; i < m; i++) lu->perm_c[i] = ra[i];
278: PetscCall(ISRestoreIndices(r, &ra));
279: }
281: /* print the control parameters */
282: if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control);
284: /* symbolic factorization of A' */
285: if (r) { /* use Petsc row ordering */
286: #if !defined(PETSC_USE_COMPLEX)
287: status = umfpack_UMF_qsymbolic(n, m, ai, aj, av, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info);
288: #else
289: status = umfpack_UMF_qsymbolic(n, m, ai, aj, NULL, NULL, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info);
290: #endif
291: } else { /* use Umfpack col ordering */
292: #if !defined(PETSC_USE_COMPLEX)
293: status = umfpack_UMF_symbolic(n, m, ai, aj, av, &lu->Symbolic, lu->Control, lu->Info);
294: #else
295: status = umfpack_UMF_symbolic(n, m, ai, aj, NULL, NULL, &lu->Symbolic, lu->Control, lu->Info);
296: #endif
297: }
298: if (status < 0) {
299: umfpack_UMF_report_info(lu->Control, lu->Info);
300: umfpack_UMF_report_status(lu->Control, status);
301: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_symbolic failed");
302: }
303: /* report sumbolic factorization of A' when Control[PRL] > 3 */
304: (void)umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control);
306: lu->flg = DIFFERENT_NONZERO_PATTERN;
307: lu->CleanUpUMFPACK = PETSC_TRUE;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: static PetscErrorCode MatView_Info_UMFPACK(Mat A, PetscViewer viewer)
312: {
313: Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;
315: PetscFunctionBegin;
316: /* check if matrix is UMFPACK type */
317: if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(PETSC_SUCCESS);
319: PetscCall(PetscViewerASCIIPrintf(viewer, "UMFPACK run parameters:\n"));
320: /* Control parameters used by reporting routiones */
321: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_PRL]: %g\n", lu->Control[UMFPACK_PRL]));
323: /* Control parameters used by symbolic factorization */
324: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_STRATEGY]: %g\n", lu->Control[UMFPACK_STRATEGY]));
325: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DENSE_COL]: %g\n", lu->Control[UMFPACK_DENSE_COL]));
326: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DENSE_ROW]: %g\n", lu->Control[UMFPACK_DENSE_ROW]));
327: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_AMD_DENSE]: %g\n", lu->Control[UMFPACK_AMD_DENSE]));
328: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_BLOCK_SIZE]: %g\n", lu->Control[UMFPACK_BLOCK_SIZE]));
329: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_FIXQ]: %g\n", lu->Control[UMFPACK_FIXQ]));
330: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_AGGRESSIVE]: %g\n", lu->Control[UMFPACK_AGGRESSIVE]));
332: /* Control parameters used by numeric factorization */
333: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_PIVOT_TOLERANCE]));
334: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]));
335: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_SCALE]: %g\n", lu->Control[UMFPACK_SCALE]));
336: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_ALLOC_INIT]: %g\n", lu->Control[UMFPACK_ALLOC_INIT]));
337: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DROPTOL]: %g\n", lu->Control[UMFPACK_DROPTOL]));
339: /* Control parameters used by solve */
340: PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_IRSTEP]: %g\n", lu->Control[UMFPACK_IRSTEP]));
342: /* mat ordering */
343: if (!lu->perm_c) PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n", UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]]));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: static PetscErrorCode MatView_UMFPACK(Mat A, PetscViewer viewer)
348: {
349: PetscBool iascii;
350: PetscViewerFormat format;
352: PetscFunctionBegin;
353: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
354: if (iascii) {
355: PetscCall(PetscViewerGetFormat(viewer, &format));
356: if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_UMFPACK(A, viewer));
357: }
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: static PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A, MatSolverType *type)
362: {
363: PetscFunctionBegin;
364: *type = MATSOLVERUMFPACK;
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*MC
369: MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers, LU, for sequential matrices
370: via the external package UMFPACK.
372: Use `./configure --download-suitesparse` to install PETSc to use UMFPACK
374: Use `-pc_type lu` `-pc_factor_mat_solver_type umfpack` to use this direct solver
376: Consult UMFPACK documentation for more information about the Control parameters
377: which correspond to the options database keys below.
379: Options Database Keys:
380: + -mat_umfpack_ordering - `CHOLMOD`, `AMD`, `GIVEN`, `METIS`, `BEST`, `NONE`
381: . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL]
382: . -mat_umfpack_strategy <AUTO> - (choose one of) `AUTO`, `UNSYMMETRIC`, `SYMMETRIC 2BY2`
383: . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
384: . -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW]
385: . -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE]
386: . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
387: . -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE]
388: . -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ]
389: . -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE]
390: . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
391: . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE]
392: . -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX
393: . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
394: . -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL]
395: - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]
397: Level: beginner
399: Note:
400: UMFPACK is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>
402: .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`
403: M*/
405: static PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A, MatFactorType ftype, Mat *F)
406: {
407: Mat B;
408: Mat_UMFPACK *lu;
409: PetscInt m = A->rmap->n, n = A->cmap->n;
411: PetscFunctionBegin;
412: /* Create the factorization matrix F */
413: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
414: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
415: PetscCall(PetscStrallocpy("umfpack", &((PetscObject)B)->type_name));
416: PetscCall(MatSetUp(B));
418: PetscCall(PetscNew(&lu));
420: B->data = lu;
421: B->ops->getinfo = MatGetInfo_External;
422: B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
423: B->ops->destroy = MatDestroy_UMFPACK;
424: B->ops->view = MatView_UMFPACK;
425: B->ops->matsolve = NULL;
427: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_umfpack));
429: B->factortype = MAT_FACTOR_LU;
430: B->assembled = PETSC_TRUE; /* required by -ksp_view */
431: B->preallocated = PETSC_TRUE;
433: PetscCall(PetscFree(B->solvertype));
434: PetscCall(PetscStrallocpy(MATSOLVERUMFPACK, &B->solvertype));
435: B->canuseordering = PETSC_TRUE;
436: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]));
438: /* initializations */
439: /* get the default control parameters */
440: umfpack_UMF_defaults(lu->Control);
441: lu->perm_c = NULL; /* use default UMFPACK col permutation */
442: lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */
444: *F = B;
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat, MatFactorType, Mat *);
449: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat, MatFactorType, Mat *);
450: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat, MatFactorType, Mat *);
451: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat, MatFactorType, Mat *);
453: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void)
454: {
455: PetscFunctionBegin;
456: PetscCall(MatSolverTypeRegister(MATSOLVERUMFPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_umfpack));
457: PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqaij_cholmod));
458: PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_cholmod));
459: PetscCall(MatSolverTypeRegister(MATSOLVERKLU, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_klu));
460: PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATSEQAIJ, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
461: if (!PetscDefined(USE_COMPLEX)) PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMAL, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
462: PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMALHERMITIAN, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }