1 2 /* 3 Provides an interface to the UMFPACK sparse solver available through SuiteSparse version 4.2.1 4 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 32bit machines [as Suitesparse_long is 32bit only] 9 10 */ 11 #include <../src/mat/impls/aij/seq/aij.h> 12 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 28 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 43 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 58 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 74 75 EXTERN_C_BEGIN 76 #include <umfpack.h> 77 EXTERN_C_END 78 79 static const char *const UmfpackOrderingTypes[] = {"CHOLMOD", "AMD", "GIVEN", "METIS", "BEST", "NONE", "USER", "UmfpackOrderingTypes", "UMFPACK_ORDERING_", 0}; 80 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; 87 88 /* Flag to clean up UMFPACK objects during Destroy */ 89 PetscBool CleanUpUMFPACK; 90 } Mat_UMFPACK; 91 92 static PetscErrorCode MatDestroy_UMFPACK(Mat A) { 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(0); 107 } 108 109 static PetscErrorCode MatSolve_UMFPACK_Private(Mat A, Vec b, Vec x, int uflag) { 110 Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data; 111 Mat_SeqAIJ *a = (Mat_SeqAIJ *)lu->A->data; 112 PetscScalar *av = a->a, *xa; 113 const PetscScalar *ba; 114 PetscInt *ai = a->i, *aj = a->j, status; 115 static PetscBool cite = PETSC_FALSE; 116 117 PetscFunctionBegin; 118 if (!A->rmap->n) PetscFunctionReturn(0); 119 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 " 120 "volume={30},\n number={2},\n pages={196--199},\n year={2004},\n publisher={ACM}\n}\n", 121 &cite)); 122 /* solve Ax = b by umfpack_*_wsolve */ 123 /* ----------------------------------*/ 124 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 } 129 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 } 142 143 PetscCall(VecRestoreArrayRead(b, &ba)); 144 PetscCall(VecRestoreArray(x, &xa)); 145 PetscFunctionReturn(0); 146 } 147 148 static PetscErrorCode MatSolve_UMFPACK(Mat A, Vec b, Vec x) { 149 PetscFunctionBegin; 150 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 151 PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_Aat)); 152 PetscFunctionReturn(0); 153 } 154 155 static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A, Vec b, Vec x) { 156 PetscFunctionBegin; 157 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 158 PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_A)); 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F, Mat A, const MatFactorInfo *info) { 163 Mat_UMFPACK *lu = (Mat_UMFPACK *)(F)->data; 164 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 165 PetscInt *ai = a->i, *aj = a->j, status; 166 PetscScalar *av = a->a; 167 168 PetscFunctionBegin; 169 if (!A->rmap->n) PetscFunctionReturn(0); 170 /* numeric factorization of A' */ 171 /* ----------------------------*/ 172 173 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) umfpack_UMF_free_numeric(&lu->Numeric); 174 #if defined(PETSC_USE_COMPLEX) 175 status = umfpack_UMF_numeric(ai, aj, (double *)av, NULL, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info); 176 #else 177 status = umfpack_UMF_numeric(ai, aj, av, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info); 178 #endif 179 if (status < 0) { 180 umfpack_UMF_report_status(lu->Control, status); 181 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_numeric failed"); 182 } 183 /* report numeric factorization of A' when Control[PRL] > 3 */ 184 (void)umfpack_UMF_report_numeric(lu->Numeric, lu->Control); 185 186 PetscCall(PetscObjectReference((PetscObject)A)); 187 PetscCall(MatDestroy(&lu->A)); 188 189 lu->A = A; 190 lu->flg = SAME_NONZERO_PATTERN; 191 lu->CleanUpUMFPACK = PETSC_TRUE; 192 F->ops->solve = MatSolve_UMFPACK; 193 F->ops->solvetranspose = MatSolveTranspose_UMFPACK; 194 PetscFunctionReturn(0); 195 } 196 197 static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) { 198 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 199 Mat_UMFPACK *lu = (Mat_UMFPACK *)(F->data); 200 PetscInt i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n, status, idx; 201 #if !defined(PETSC_USE_COMPLEX) 202 PetscScalar *av = a->a; 203 #endif 204 const PetscInt *ra; 205 const char *strategy[] = {"AUTO", "UNSYMMETRIC", "SYMMETRIC"}; 206 const char *scale[] = {"NONE", "SUM", "MAX"}; 207 PetscBool flg; 208 209 PetscFunctionBegin; 210 (F)->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK; 211 if (!n) PetscFunctionReturn(0); 212 213 /* Set options to F */ 214 PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "UMFPACK Options", "Mat"); 215 /* Control parameters used by reporting routiones */ 216 PetscCall(PetscOptionsReal("-mat_umfpack_prl", "Control[UMFPACK_PRL]", "None", lu->Control[UMFPACK_PRL], &lu->Control[UMFPACK_PRL], NULL)); 217 218 /* Control parameters for symbolic factorization */ 219 PetscCall(PetscOptionsEList("-mat_umfpack_strategy", "ordering and pivoting strategy", "None", strategy, 3, strategy[0], &idx, &flg)); 220 if (flg) { 221 switch (idx) { 222 case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break; 223 case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break; 224 case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break; 225 } 226 } 227 PetscCall(PetscOptionsEList("-mat_umfpack_ordering", "Internal ordering method", "None", UmfpackOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(UmfpackOrderingTypes), UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]], &idx, &flg)); 228 if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx; 229 PetscCall(PetscOptionsReal("-mat_umfpack_dense_col", "Control[UMFPACK_DENSE_COL]", "None", lu->Control[UMFPACK_DENSE_COL], &lu->Control[UMFPACK_DENSE_COL], NULL)); 230 PetscCall(PetscOptionsReal("-mat_umfpack_dense_row", "Control[UMFPACK_DENSE_ROW]", "None", lu->Control[UMFPACK_DENSE_ROW], &lu->Control[UMFPACK_DENSE_ROW], NULL)); 231 PetscCall(PetscOptionsReal("-mat_umfpack_amd_dense", "Control[UMFPACK_AMD_DENSE]", "None", lu->Control[UMFPACK_AMD_DENSE], &lu->Control[UMFPACK_AMD_DENSE], NULL)); 232 PetscCall(PetscOptionsReal("-mat_umfpack_block_size", "Control[UMFPACK_BLOCK_SIZE]", "None", lu->Control[UMFPACK_BLOCK_SIZE], &lu->Control[UMFPACK_BLOCK_SIZE], NULL)); 233 PetscCall(PetscOptionsReal("-mat_umfpack_fixq", "Control[UMFPACK_FIXQ]", "None", lu->Control[UMFPACK_FIXQ], &lu->Control[UMFPACK_FIXQ], NULL)); 234 PetscCall(PetscOptionsReal("-mat_umfpack_aggressive", "Control[UMFPACK_AGGRESSIVE]", "None", lu->Control[UMFPACK_AGGRESSIVE], &lu->Control[UMFPACK_AGGRESSIVE], NULL)); 235 236 /* Control parameters used by numeric factorization */ 237 PetscCall(PetscOptionsReal("-mat_umfpack_pivot_tolerance", "Control[UMFPACK_PIVOT_TOLERANCE]", "None", lu->Control[UMFPACK_PIVOT_TOLERANCE], &lu->Control[UMFPACK_PIVOT_TOLERANCE], NULL)); 238 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)); 239 PetscCall(PetscOptionsEList("-mat_umfpack_scale", "Control[UMFPACK_SCALE]", "None", scale, 3, scale[0], &idx, &flg)); 240 if (flg) { 241 switch (idx) { 242 case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break; 243 case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break; 244 case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break; 245 } 246 } 247 PetscCall(PetscOptionsReal("-mat_umfpack_alloc_init", "Control[UMFPACK_ALLOC_INIT]", "None", lu->Control[UMFPACK_ALLOC_INIT], &lu->Control[UMFPACK_ALLOC_INIT], NULL)); 248 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)); 249 PetscCall(PetscOptionsReal("-mat_umfpack_droptol", "Control[UMFPACK_DROPTOL]", "None", lu->Control[UMFPACK_DROPTOL], &lu->Control[UMFPACK_DROPTOL], NULL)); 250 251 /* Control parameters used by solve */ 252 PetscCall(PetscOptionsReal("-mat_umfpack_irstep", "Control[UMFPACK_IRSTEP]", "None", lu->Control[UMFPACK_IRSTEP], &lu->Control[UMFPACK_IRSTEP], NULL)); 253 PetscOptionsEnd(); 254 255 if (r) { 256 PetscCall(ISGetIndices(r, &ra)); 257 PetscCall(PetscMalloc1(m, &lu->perm_c)); 258 /* we cannot simply memcpy on 64 bit archs */ 259 for (i = 0; i < m; i++) lu->perm_c[i] = ra[i]; 260 PetscCall(ISRestoreIndices(r, &ra)); 261 } 262 263 /* print the control parameters */ 264 if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control); 265 266 /* symbolic factorization of A' */ 267 /* ---------------------------------------------------------------------- */ 268 if (r) { /* use Petsc row ordering */ 269 #if !defined(PETSC_USE_COMPLEX) 270 status = umfpack_UMF_qsymbolic(n, m, ai, aj, av, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info); 271 #else 272 status = umfpack_UMF_qsymbolic(n, m, ai, aj, NULL, NULL, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info); 273 #endif 274 } else { /* use Umfpack col ordering */ 275 #if !defined(PETSC_USE_COMPLEX) 276 status = umfpack_UMF_symbolic(n, m, ai, aj, av, &lu->Symbolic, lu->Control, lu->Info); 277 #else 278 status = umfpack_UMF_symbolic(n, m, ai, aj, NULL, NULL, &lu->Symbolic, lu->Control, lu->Info); 279 #endif 280 } 281 if (status < 0) { 282 umfpack_UMF_report_info(lu->Control, lu->Info); 283 umfpack_UMF_report_status(lu->Control, status); 284 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_symbolic failed"); 285 } 286 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 287 (void)umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control); 288 289 lu->flg = DIFFERENT_NONZERO_PATTERN; 290 lu->CleanUpUMFPACK = PETSC_TRUE; 291 PetscFunctionReturn(0); 292 } 293 294 static PetscErrorCode MatView_Info_UMFPACK(Mat A, PetscViewer viewer) { 295 Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data; 296 297 PetscFunctionBegin; 298 /* check if matrix is UMFPACK type */ 299 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0); 300 301 PetscCall(PetscViewerASCIIPrintf(viewer, "UMFPACK run parameters:\n")); 302 /* Control parameters used by reporting routiones */ 303 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_PRL]: %g\n", lu->Control[UMFPACK_PRL])); 304 305 /* Control parameters used by symbolic factorization */ 306 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_STRATEGY]: %g\n", lu->Control[UMFPACK_STRATEGY])); 307 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DENSE_COL]: %g\n", lu->Control[UMFPACK_DENSE_COL])); 308 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DENSE_ROW]: %g\n", lu->Control[UMFPACK_DENSE_ROW])); 309 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_AMD_DENSE]: %g\n", lu->Control[UMFPACK_AMD_DENSE])); 310 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_BLOCK_SIZE]: %g\n", lu->Control[UMFPACK_BLOCK_SIZE])); 311 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_FIXQ]: %g\n", lu->Control[UMFPACK_FIXQ])); 312 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_AGGRESSIVE]: %g\n", lu->Control[UMFPACK_AGGRESSIVE])); 313 314 /* Control parameters used by numeric factorization */ 315 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_PIVOT_TOLERANCE])); 316 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE])); 317 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_SCALE]: %g\n", lu->Control[UMFPACK_SCALE])); 318 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_ALLOC_INIT]: %g\n", lu->Control[UMFPACK_ALLOC_INIT])); 319 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DROPTOL]: %g\n", lu->Control[UMFPACK_DROPTOL])); 320 321 /* Control parameters used by solve */ 322 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_IRSTEP]: %g\n", lu->Control[UMFPACK_IRSTEP])); 323 324 /* mat ordering */ 325 if (!lu->perm_c) PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n", UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]])); 326 PetscFunctionReturn(0); 327 } 328 329 static PetscErrorCode MatView_UMFPACK(Mat A, PetscViewer viewer) { 330 PetscBool iascii; 331 PetscViewerFormat format; 332 333 PetscFunctionBegin; 334 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 335 if (iascii) { 336 PetscCall(PetscViewerGetFormat(viewer, &format)); 337 if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_UMFPACK(A, viewer)); 338 } 339 PetscFunctionReturn(0); 340 } 341 342 PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A, MatSolverType *type) { 343 PetscFunctionBegin; 344 *type = MATSOLVERUMFPACK; 345 PetscFunctionReturn(0); 346 } 347 348 /*MC 349 MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers, LU, for sequential matrices 350 via the external package UMFPACK. 351 352 Use ./configure --download-suitesparse to install PETSc to use UMFPACK 353 354 Use -pc_type lu -pc_factor_mat_solver_type umfpack to use this direct solver 355 356 Consult UMFPACK documentation for more information about the Control parameters 357 which correspond to the options database keys below. 358 359 Options Database Keys: 360 + -mat_umfpack_ordering - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE 361 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 362 . -mat_umfpack_strategy <AUTO> - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2 363 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 364 . -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW] 365 . -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE] 366 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 367 . -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE] 368 . -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ] 369 . -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE] 370 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 371 . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE] 372 . -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX 373 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 374 . -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL] 375 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 376 377 Level: beginner 378 379 Note: UMFPACK is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 380 381 .seealso: `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType` 382 M*/ 383 384 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A, MatFactorType ftype, Mat *F) { 385 Mat B; 386 Mat_UMFPACK *lu; 387 PetscInt m = A->rmap->n, n = A->cmap->n; 388 389 PetscFunctionBegin; 390 /* Create the factorization matrix F */ 391 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 392 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n)); 393 PetscCall(PetscStrallocpy("umfpack", &((PetscObject)B)->type_name)); 394 PetscCall(MatSetUp(B)); 395 396 PetscCall(PetscNewLog(B, &lu)); 397 398 B->data = lu; 399 B->ops->getinfo = MatGetInfo_External; 400 B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK; 401 B->ops->destroy = MatDestroy_UMFPACK; 402 B->ops->view = MatView_UMFPACK; 403 B->ops->matsolve = NULL; 404 405 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_umfpack)); 406 407 B->factortype = MAT_FACTOR_LU; 408 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 409 B->preallocated = PETSC_TRUE; 410 411 PetscCall(PetscFree(B->solvertype)); 412 PetscCall(PetscStrallocpy(MATSOLVERUMFPACK, &B->solvertype)); 413 B->canuseordering = PETSC_TRUE; 414 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU])); 415 416 /* initializations */ 417 /* ------------------------------------------------*/ 418 /* get the default control parameters */ 419 umfpack_UMF_defaults(lu->Control); 420 lu->perm_c = NULL; /* use defaul UMFPACK col permutation */ 421 lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */ 422 423 *F = B; 424 PetscFunctionReturn(0); 425 } 426 427 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat, MatFactorType, Mat *); 428 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat, MatFactorType, Mat *); 429 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat, MatFactorType, Mat *); 430 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat, MatFactorType, Mat *); 431 432 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void) { 433 PetscFunctionBegin; 434 PetscCall(MatSolverTypeRegister(MATSOLVERUMFPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_umfpack)); 435 PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqaij_cholmod)); 436 PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_cholmod)); 437 PetscCall(MatSolverTypeRegister(MATSOLVERKLU, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_klu)); 438 PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATSEQAIJ, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr)); 439 if (!PetscDefined(USE_COMPLEX)) PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMAL, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr)); 440 PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMALHERMITIAN, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr)); 441 PetscFunctionReturn(0); 442 } 443