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