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