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 32-bit pointer systems [as Suitesparse_long is 32-bit 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:   PetscFunctionBegin;
 97:   if (lu->CleanUpUMFPACK) {
 98:     umfpack_UMF_free_symbolic(&lu->Symbolic);
 99:     umfpack_UMF_free_numeric(&lu->Numeric);
100:     PetscCall(PetscFree(lu->Wi));
101:     PetscCall(PetscFree(lu->W));
102:     PetscCall(PetscFree(lu->perm_c));
103:   }
104:   PetscCall(MatDestroy(&lu->A));
105:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
106:   PetscCall(PetscFree(A->data));
107:   PetscFunctionReturn(PETSC_SUCCESS);
108: }

110: static PetscErrorCode MatSolve_UMFPACK_Private(Mat A, Vec b, Vec x, int uflag)
111: {
112:   Mat_UMFPACK       *lu = (Mat_UMFPACK *)A->data;
113:   Mat_SeqAIJ        *a  = (Mat_SeqAIJ *)lu->A->data;
114:   PetscScalar       *av = a->a, *xa;
115:   const PetscScalar *ba;
116:   PetscInt          *ai = a->i, *aj = a->j, 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, status;
170:   PetscScalar *av = a->a;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

395:    Level: beginner

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

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

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

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

416:   PetscCall(PetscNew(&lu));

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

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

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

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

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

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

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

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