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 { 94 Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data; 95 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 } 109 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; 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 127 if (!lu->Wi) { /* first time, allocate working space for wsolve */ 128 PetscCall(PetscMalloc1(A->rmap->n, &lu->Wi)); 129 PetscCall(PetscMalloc1(5 * A->rmap->n, &lu->W)); 130 } 131 132 PetscCall(VecGetArrayRead(b, &ba)); 133 PetscCall(VecGetArray(x, &xa)); 134 #if defined(PETSC_USE_COMPLEX) 135 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); 136 #else 137 status = umfpack_UMF_wsolve(uflag, ai, aj, av, xa, ba, lu->Numeric, lu->Control, lu->Info, lu->Wi, lu->W); 138 #endif 139 umfpack_UMF_report_info(lu->Control, lu->Info); 140 if (status < 0) { 141 umfpack_UMF_report_status(lu->Control, status); 142 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_wsolve failed"); 143 } 144 145 PetscCall(VecRestoreArrayRead(b, &ba)); 146 PetscCall(VecRestoreArray(x, &xa)); 147 PetscFunctionReturn(PETSC_SUCCESS); 148 } 149 150 static PetscErrorCode MatSolve_UMFPACK(Mat A, Vec b, Vec x) 151 { 152 PetscFunctionBegin; 153 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 154 PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_Aat)); 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A, Vec b, Vec x) 159 { 160 PetscFunctionBegin; 161 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 162 PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_A)); 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F, Mat A, const MatFactorInfo *info) 167 { 168 Mat_UMFPACK *lu = (Mat_UMFPACK *)(F)->data; 169 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 170 PetscInt *ai = a->i, *aj = a->j, status; 171 PetscScalar *av = a->a; 172 173 PetscFunctionBegin; 174 if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS); 175 /* numeric factorization of A' */ 176 /* ----------------------------*/ 177 178 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) umfpack_UMF_free_numeric(&lu->Numeric); 179 #if defined(PETSC_USE_COMPLEX) 180 status = umfpack_UMF_numeric(ai, aj, (double *)av, NULL, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info); 181 #else 182 status = umfpack_UMF_numeric(ai, aj, av, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info); 183 #endif 184 if (status < 0) { 185 umfpack_UMF_report_status(lu->Control, status); 186 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_numeric failed"); 187 } 188 /* report numeric factorization of A' when Control[PRL] > 3 */ 189 (void)umfpack_UMF_report_numeric(lu->Numeric, lu->Control); 190 191 PetscCall(PetscObjectReference((PetscObject)A)); 192 PetscCall(MatDestroy(&lu->A)); 193 194 lu->A = A; 195 lu->flg = SAME_NONZERO_PATTERN; 196 lu->CleanUpUMFPACK = PETSC_TRUE; 197 F->ops->solve = MatSolve_UMFPACK; 198 F->ops->solvetranspose = MatSolveTranspose_UMFPACK; 199 PetscFunctionReturn(PETSC_SUCCESS); 200 } 201 202 static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) 203 { 204 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 205 Mat_UMFPACK *lu = (Mat_UMFPACK *)(F->data); 206 PetscInt i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n, status, idx; 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 /* ---------------------------------------------------------------------- */ 286 if (r) { /* use Petsc row ordering */ 287 #if !defined(PETSC_USE_COMPLEX) 288 status = umfpack_UMF_qsymbolic(n, m, ai, aj, av, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info); 289 #else 290 status = umfpack_UMF_qsymbolic(n, m, ai, aj, NULL, NULL, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info); 291 #endif 292 } else { /* use Umfpack col ordering */ 293 #if !defined(PETSC_USE_COMPLEX) 294 status = umfpack_UMF_symbolic(n, m, ai, aj, av, &lu->Symbolic, lu->Control, lu->Info); 295 #else 296 status = umfpack_UMF_symbolic(n, m, ai, aj, NULL, NULL, &lu->Symbolic, lu->Control, lu->Info); 297 #endif 298 } 299 if (status < 0) { 300 umfpack_UMF_report_info(lu->Control, lu->Info); 301 umfpack_UMF_report_status(lu->Control, status); 302 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_symbolic failed"); 303 } 304 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 305 (void)umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control); 306 307 lu->flg = DIFFERENT_NONZERO_PATTERN; 308 lu->CleanUpUMFPACK = PETSC_TRUE; 309 PetscFunctionReturn(PETSC_SUCCESS); 310 } 311 312 static PetscErrorCode MatView_Info_UMFPACK(Mat A, PetscViewer viewer) 313 { 314 Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data; 315 316 PetscFunctionBegin; 317 /* check if matrix is UMFPACK type */ 318 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(PETSC_SUCCESS); 319 320 PetscCall(PetscViewerASCIIPrintf(viewer, "UMFPACK run parameters:\n")); 321 /* Control parameters used by reporting routiones */ 322 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_PRL]: %g\n", lu->Control[UMFPACK_PRL])); 323 324 /* Control parameters used by symbolic factorization */ 325 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_STRATEGY]: %g\n", lu->Control[UMFPACK_STRATEGY])); 326 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DENSE_COL]: %g\n", lu->Control[UMFPACK_DENSE_COL])); 327 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DENSE_ROW]: %g\n", lu->Control[UMFPACK_DENSE_ROW])); 328 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_AMD_DENSE]: %g\n", lu->Control[UMFPACK_AMD_DENSE])); 329 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_BLOCK_SIZE]: %g\n", lu->Control[UMFPACK_BLOCK_SIZE])); 330 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_FIXQ]: %g\n", lu->Control[UMFPACK_FIXQ])); 331 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_AGGRESSIVE]: %g\n", lu->Control[UMFPACK_AGGRESSIVE])); 332 333 /* Control parameters used by numeric factorization */ 334 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_PIVOT_TOLERANCE])); 335 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE])); 336 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_SCALE]: %g\n", lu->Control[UMFPACK_SCALE])); 337 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_ALLOC_INIT]: %g\n", lu->Control[UMFPACK_ALLOC_INIT])); 338 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_DROPTOL]: %g\n", lu->Control[UMFPACK_DROPTOL])); 339 340 /* Control parameters used by solve */ 341 PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_IRSTEP]: %g\n", lu->Control[UMFPACK_IRSTEP])); 342 343 /* mat ordering */ 344 if (!lu->perm_c) PetscCall(PetscViewerASCIIPrintf(viewer, " Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n", UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]])); 345 PetscFunctionReturn(PETSC_SUCCESS); 346 } 347 348 static PetscErrorCode MatView_UMFPACK(Mat A, PetscViewer viewer) 349 { 350 PetscBool iascii; 351 PetscViewerFormat format; 352 353 PetscFunctionBegin; 354 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 355 if (iascii) { 356 PetscCall(PetscViewerGetFormat(viewer, &format)); 357 if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_UMFPACK(A, viewer)); 358 } 359 PetscFunctionReturn(PETSC_SUCCESS); 360 } 361 362 PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A, MatSolverType *type) 363 { 364 PetscFunctionBegin; 365 *type = MATSOLVERUMFPACK; 366 PetscFunctionReturn(PETSC_SUCCESS); 367 } 368 369 /*MC 370 MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers, LU, for sequential matrices 371 via the external package UMFPACK. 372 373 Use ./configure --download-suitesparse to install PETSc to use UMFPACK 374 375 Use -pc_type lu -pc_factor_mat_solver_type umfpack to use this direct solver 376 377 Consult UMFPACK documentation for more information about the Control parameters 378 which correspond to the options database keys below. 379 380 Options Database Keys: 381 + -mat_umfpack_ordering - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE 382 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 383 . -mat_umfpack_strategy <AUTO> - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2 384 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 385 . -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW] 386 . -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE] 387 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 388 . -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE] 389 . -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ] 390 . -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE] 391 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 392 . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE] 393 . -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX 394 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 395 . -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL] 396 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 397 398 Level: beginner 399 400 Note: UMFPACK is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 401 402 .seealso: `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType` 403 M*/ 404 405 PETSC_EXTERN 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 /* ------------------------------------------------*/ 440 /* get the default control parameters */ 441 umfpack_UMF_defaults(lu->Control); 442 lu->perm_c = NULL; /* use defaul UMFPACK col permutation */ 443 lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */ 444 445 *F = B; 446 PetscFunctionReturn(PETSC_SUCCESS); 447 } 448 449 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat, MatFactorType, Mat *); 450 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat, MatFactorType, Mat *); 451 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat, MatFactorType, Mat *); 452 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat, MatFactorType, Mat *); 453 454 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void) 455 { 456 PetscFunctionBegin; 457 PetscCall(MatSolverTypeRegister(MATSOLVERUMFPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_umfpack)); 458 PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqaij_cholmod)); 459 PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_cholmod)); 460 PetscCall(MatSolverTypeRegister(MATSOLVERKLU, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_klu)); 461 PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATSEQAIJ, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr)); 462 if (!PetscDefined(USE_COMPLEX)) PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMAL, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr)); 463 PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMALHERMITIAN, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr)); 464 PetscFunctionReturn(PETSC_SUCCESS); 465 } 466