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