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