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(PetscFree(A->data)); 106 PetscFunctionReturn(0); 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,status; 116 static PetscBool cite = PETSC_FALSE; 117 118 PetscFunctionBegin; 119 if (!A->rmap->n) PetscFunctionReturn(0); 120 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)); 121 /* solve Ax = b by umfpack_*_wsolve */ 122 /* ----------------------------------*/ 123 124 if (!lu->Wi) { /* first time, allocate working space for wsolve */ 125 PetscCall(PetscMalloc1(A->rmap->n,&lu->Wi)); 126 PetscCall(PetscMalloc1(5*A->rmap->n,&lu->W)); 127 } 128 129 PetscCall(VecGetArrayRead(b,&ba)); 130 PetscCall(VecGetArray(x,&xa)); 131 #if defined(PETSC_USE_COMPLEX) 132 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); 133 #else 134 status = umfpack_UMF_wsolve(uflag,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 135 #endif 136 umfpack_UMF_report_info(lu->Control, lu->Info); 137 if (status < 0) { 138 umfpack_UMF_report_status(lu->Control, status); 139 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_wsolve failed"); 140 } 141 142 PetscCall(VecRestoreArrayRead(b,&ba)); 143 PetscCall(VecRestoreArray(x,&xa)); 144 PetscFunctionReturn(0); 145 } 146 147 static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x) 148 { 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 { 157 PetscFunctionBegin; 158 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 159 PetscCall(MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A)); 160 PetscFunctionReturn(0); 161 } 162 163 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info) 164 { 165 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F)->data; 166 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 167 PetscInt *ai = a->i,*aj=a->j,status; 168 PetscScalar *av = a->a; 169 170 PetscFunctionBegin; 171 if (!A->rmap->n) PetscFunctionReturn(0); 172 /* numeric factorization of A' */ 173 /* ----------------------------*/ 174 175 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) { 176 umfpack_UMF_free_numeric(&lu->Numeric); 177 } 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(0); 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; 206 #if !defined(PETSC_USE_COMPLEX) 207 PetscScalar *av = a->a; 208 #endif 209 const PetscInt *ra; 210 PetscInt status; 211 212 PetscFunctionBegin; 213 (F)->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK; 214 if (!n) PetscFunctionReturn(0); 215 if (r) { 216 PetscCall(ISGetIndices(r,&ra)); 217 PetscCall(PetscMalloc1(m,&lu->perm_c)); 218 /* we cannot simply memcpy on 64 bit archs */ 219 for (i = 0; i < m; i++) lu->perm_c[i] = ra[i]; 220 PetscCall(ISRestoreIndices(r,&ra)); 221 } 222 223 /* print the control parameters */ 224 if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control); 225 226 /* symbolic factorization of A' */ 227 /* ---------------------------------------------------------------------- */ 228 if (r) { /* use Petsc row ordering */ 229 #if !defined(PETSC_USE_COMPLEX) 230 status = umfpack_UMF_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 231 #else 232 status = umfpack_UMF_qsymbolic(n,m,ai,aj,NULL,NULL,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 233 #endif 234 } else { /* use Umfpack col ordering */ 235 #if !defined(PETSC_USE_COMPLEX) 236 status = umfpack_UMF_symbolic(n,m,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info); 237 #else 238 status = umfpack_UMF_symbolic(n,m,ai,aj,NULL,NULL,&lu->Symbolic,lu->Control,lu->Info); 239 #endif 240 } 241 if (status < 0) { 242 umfpack_UMF_report_info(lu->Control, lu->Info); 243 umfpack_UMF_report_status(lu->Control, status); 244 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_symbolic failed"); 245 } 246 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 247 (void) umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control); 248 249 lu->flg = DIFFERENT_NONZERO_PATTERN; 250 lu->CleanUpUMFPACK = PETSC_TRUE; 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode MatView_Info_UMFPACK(Mat A,PetscViewer viewer) 255 { 256 Mat_UMFPACK *lu= (Mat_UMFPACK*)A->data; 257 258 PetscFunctionBegin; 259 /* check if matrix is UMFPACK type */ 260 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0); 261 262 PetscCall(PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n")); 263 /* Control parameters used by reporting routiones */ 264 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL])); 265 266 /* Control parameters used by symbolic factorization */ 267 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY])); 268 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL])); 269 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW])); 270 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE])); 271 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE])); 272 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ])); 273 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE])); 274 275 /* Control parameters used by numeric factorization */ 276 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE])); 277 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE])); 278 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE])); 279 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT])); 280 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL])); 281 282 /* Control parameters used by solve */ 283 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP])); 284 285 /* mat ordering */ 286 if (!lu->perm_c) { 287 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n",UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]])); 288 } 289 PetscFunctionReturn(0); 290 } 291 292 static PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer) 293 { 294 PetscBool iascii; 295 PetscViewerFormat format; 296 297 PetscFunctionBegin; 298 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 299 if (iascii) { 300 PetscCall(PetscViewerGetFormat(viewer,&format)); 301 if (format == PETSC_VIEWER_ASCII_INFO) { 302 PetscCall(MatView_Info_UMFPACK(A,viewer)); 303 } 304 } 305 PetscFunctionReturn(0); 306 } 307 308 PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A,MatSolverType *type) 309 { 310 PetscFunctionBegin; 311 *type = MATSOLVERUMFPACK; 312 PetscFunctionReturn(0); 313 } 314 315 /*MC 316 MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices 317 via the external package UMFPACK. 318 319 Use ./configure --download-suitesparse to install PETSc to use UMFPACK 320 321 Use -pc_type lu -pc_factor_mat_solver_type umfpack to use this direct solver 322 323 Consult UMFPACK documentation for more information about the Control parameters 324 which correspond to the options database keys below. 325 326 Options Database Keys: 327 + -mat_umfpack_ordering - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE 328 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 329 . -mat_umfpack_strategy <AUTO> - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2 330 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 331 . -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW] 332 . -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE] 333 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 334 . -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE] 335 . -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ] 336 . -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE] 337 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 338 . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE] 339 . -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX 340 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 341 . -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL] 342 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 343 344 Level: beginner 345 346 Note: UMFPACK is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 347 348 .seealso: PCLU, MATSOLVERSUPERLU, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType 349 M*/ 350 351 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F) 352 { 353 Mat B; 354 Mat_UMFPACK *lu; 355 PetscErrorCode ierr; 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 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"UMFPACK Options","Mat");PetscCall(ierr); 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,sizeof(UmfpackOrderingTypes)/sizeof(UmfpackOrderingTypes[0]),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 ierr = PetscOptionsEnd();PetscCall(ierr); 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