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