1 /*
2 Provides an interface to the UMFPACK sparse solver available through SuiteSparse version 4.2.1
3
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]
8
9 */
10 #include <../src/mat/impls/aij/seq/aij.h>
11
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
27
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
42
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
57
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
73
74 EXTERN_C_BEGIN
75 #include <umfpack.h>
76 EXTERN_C_END
77
78 static const char *const UmfpackOrderingTypes[] = {"CHOLMOD", "AMD", "GIVEN", "METIS", "BEST", "NONE", "USER", "UmfpackOrderingTypes", "UMFPACK_ORDERING_", 0};
79
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;
86
87 /* Flag to clean up UMFPACK objects during Destroy */
88 PetscBool CleanUpUMFPACK;
89 } Mat_UMFPACK;
90
MatDestroy_UMFPACK(Mat A)91 static PetscErrorCode MatDestroy_UMFPACK(Mat A)
92 {
93 Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;
94
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 }
108
MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag)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;
118
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 */
125
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 }
130
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 }
143
144 PetscCall(VecRestoreArrayRead(b, &ba));
145 PetscCall(VecRestoreArray(x, &xa));
146 PetscFunctionReturn(PETSC_SUCCESS);
147 }
148
MatSolve_UMFPACK(Mat A,Vec b,Vec x)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 }
156
MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x)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 }
164
MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo * info)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;
172
173 PetscFunctionBegin;
174 if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS);
175 /* numeric factorization of A' */
176
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);
189
190 PetscCall(PetscObjectReference((PetscObject)A));
191 PetscCall(MatDestroy(&lu->A));
192
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 }
200
MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)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;
214
215 PetscFunctionBegin;
216 F->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK;
217 if (!n) PetscFunctionReturn(PETSC_SUCCESS);
218
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));
223
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));
247
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));
268
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();
272
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 }
280
281 /* print the control parameters */
282 if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control);
283
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);
305
306 lu->flg = DIFFERENT_NONZERO_PATTERN;
307 lu->CleanUpUMFPACK = PETSC_TRUE;
308 PetscFunctionReturn(PETSC_SUCCESS);
309 }
310
MatView_Info_UMFPACK(Mat A,PetscViewer viewer)311 static PetscErrorCode MatView_Info_UMFPACK(Mat A, PetscViewer viewer)
312 {
313 Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;
314
315 PetscFunctionBegin;
316 /* check if matrix is UMFPACK type */
317 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(PETSC_SUCCESS);
318
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]));
322
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]));
331
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]));
338
339 /* Control parameters used by solve */
340 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_IRSTEP]: %g\n", lu->Control[UMFPACK_IRSTEP]));
341
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 }
346
MatView_UMFPACK(Mat A,PetscViewer viewer)347 static PetscErrorCode MatView_UMFPACK(Mat A, PetscViewer viewer)
348 {
349 PetscBool isascii;
350 PetscViewerFormat format;
351
352 PetscFunctionBegin;
353 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
354 if (isascii) {
355 PetscCall(PetscViewerGetFormat(viewer, &format));
356 if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_UMFPACK(A, viewer));
357 }
358 PetscFunctionReturn(PETSC_SUCCESS);
359 }
360
MatFactorGetSolverType_seqaij_umfpack(Mat A,MatSolverType * type)361 static PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A, MatSolverType *type)
362 {
363 PetscFunctionBegin;
364 *type = MATSOLVERUMFPACK;
365 PetscFunctionReturn(PETSC_SUCCESS);
366 }
367
368 /*MC
369 MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers, LU, for sequential matrices
370 via the external package UMFPACK.
371
372 Use `./configure --download-suitesparse` to install PETSc to use UMFPACK
373
374 Use `-pc_type lu` `-pc_factor_mat_solver_type umfpack` to use this direct solver
375
376 Consult UMFPACK documentation for more information about the Control parameters
377 which correspond to the options database keys below.
378
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]
396
397 Level: beginner
398
399 Note:
400 UMFPACK is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>
401
402 .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`
403 M*/
404
MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat * F)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;
410
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));
417
418 PetscCall(PetscNew(&lu));
419
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;
426
427 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_umfpack));
428
429 B->factortype = MAT_FACTOR_LU;
430 B->assembled = PETSC_TRUE; /* required by -ksp_view */
431 B->preallocated = PETSC_TRUE;
432
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]));
437
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 */
443
444 *F = B;
445 PetscFunctionReturn(PETSC_SUCCESS);
446 }
447
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 *);
452
MatSolverTypeRegister_SuiteSparse(void)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 }
465