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, status;
116:   static PetscBool   cite = PETSC_FALSE;

118:   PetscFunctionBegin;
119:   if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS);
120:   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  "
121:                                    "volume={30},\n  number={2},\n  pages={196--199},\n  year={2004},\n  publisher={ACM}\n}\n",
122:                                    &cite));
123:   /* solve Ax = b by umfpack_*_wsolve */

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

130:   PetscCall(VecGetArrayRead(b, &ba));
131:   PetscCall(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:   PetscCall(VecRestoreArrayRead(b, &ba));
144:   PetscCall(VecRestoreArray(x, &xa));
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

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

156: static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A, Vec b, Vec x)
157: {
158:   PetscFunctionBegin;
159:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
160:   PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_A));
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

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

171:   PetscFunctionBegin;
172:   if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS);
173:   /* numeric factorization of A' */

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

188:   PetscCall(PetscObjectReference((PetscObject)A));
189:   PetscCall(MatDestroy(&lu->A));

191:   lu->A                  = A;
192:   lu->flg                = SAME_NONZERO_PATTERN;
193:   lu->CleanUpUMFPACK     = PETSC_TRUE;
194:   F->ops->solve          = MatSolve_UMFPACK;
195:   F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

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

212:   PetscFunctionBegin;
213:   (F)->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK;
214:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

266:   /* Control parameters used by solve */
267:   PetscCall(PetscOptionsReal("-mat_umfpack_irstep", "Control[UMFPACK_IRSTEP]", "None", lu->Control[UMFPACK_IRSTEP], &lu->Control[UMFPACK_IRSTEP], NULL));
268:   PetscOptionsEnd();

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

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

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

303:   lu->flg            = DIFFERENT_NONZERO_PATTERN;
304:   lu->CleanUpUMFPACK = PETSC_TRUE;
305:   PetscFunctionReturn(PETSC_SUCCESS);
306: }

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

312:   PetscFunctionBegin;
313:   /* check if matrix is UMFPACK type */
314:   if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

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

339:   /* mat ordering */
340:   if (!lu->perm_c) PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n", UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]]));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }

344: static PetscErrorCode MatView_UMFPACK(Mat A, PetscViewer viewer)
345: {
346:   PetscBool         iascii;
347:   PetscViewerFormat format;

349:   PetscFunctionBegin;
350:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
351:   if (iascii) {
352:     PetscCall(PetscViewerGetFormat(viewer, &format));
353:     if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_UMFPACK(A, viewer));
354:   }
355:   PetscFunctionReturn(PETSC_SUCCESS);
356: }

358: static PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A, MatSolverType *type)
359: {
360:   PetscFunctionBegin;
361:   *type = MATSOLVERUMFPACK;
362:   PetscFunctionReturn(PETSC_SUCCESS);
363: }

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

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

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

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

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

394:    Level: beginner

396:    Note:
397:    UMFPACK is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>

399: .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`
400: M*/

402: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A, MatFactorType ftype, Mat *F)
403: {
404:   Mat          B;
405:   Mat_UMFPACK *lu;
406:   PetscInt     m = A->rmap->n, n = A->cmap->n;

408:   PetscFunctionBegin;
409:   /* Create the factorization matrix F */
410:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
411:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
412:   PetscCall(PetscStrallocpy("umfpack", &((PetscObject)B)->type_name));
413:   PetscCall(MatSetUp(B));

415:   PetscCall(PetscNew(&lu));

417:   B->data                  = lu;
418:   B->ops->getinfo          = MatGetInfo_External;
419:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
420:   B->ops->destroy          = MatDestroy_UMFPACK;
421:   B->ops->view             = MatView_UMFPACK;
422:   B->ops->matsolve         = NULL;

424:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_umfpack));

426:   B->factortype   = MAT_FACTOR_LU;
427:   B->assembled    = PETSC_TRUE; /* required by -ksp_view */
428:   B->preallocated = PETSC_TRUE;

430:   PetscCall(PetscFree(B->solvertype));
431:   PetscCall(PetscStrallocpy(MATSOLVERUMFPACK, &B->solvertype));
432:   B->canuseordering = PETSC_TRUE;
433:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]));

435:   /* initializations */
436:   /* get the default control parameters */
437:   umfpack_UMF_defaults(lu->Control);
438:   lu->perm_c                  = NULL; /* use default UMFPACK col permutation */
439:   lu->Control[UMFPACK_IRSTEP] = 0;    /* max num of iterative refinement steps to attempt */

441:   *F = B;
442:   PetscFunctionReturn(PETSC_SUCCESS);
443: }

445: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat, MatFactorType, Mat *);
446: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat, MatFactorType, Mat *);
447: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat, MatFactorType, Mat *);
448: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat, MatFactorType, Mat *);

450: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void)
451: {
452:   PetscFunctionBegin;
453:   PetscCall(MatSolverTypeRegister(MATSOLVERUMFPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_umfpack));
454:   PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqaij_cholmod));
455:   PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_cholmod));
456:   PetscCall(MatSolverTypeRegister(MATSOLVERKLU, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_klu));
457:   PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATSEQAIJ, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
458:   if (!PetscDefined(USE_COMPLEX)) PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMAL, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
459:   PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMALHERMITIAN, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }