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