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 UF_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 UMFPACK UL_Long version MUST be built with 64 bit integers when used. 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 #define umfpack_UMF_wsolve umfpack_zl_wsolve 18 #define umfpack_UMF_numeric umfpack_zl_numeric 19 #define umfpack_UMF_report_numeric umfpack_zl_report_numeric 20 #define umfpack_UMF_report_control umfpack_zl_report_control 21 #define umfpack_UMF_report_status umfpack_zl_report_status 22 #define umfpack_UMF_report_info umfpack_zl_report_info 23 #define umfpack_UMF_report_symbolic umfpack_zl_report_symbolic 24 #define umfpack_UMF_qsymbolic umfpack_zl_qsymbolic 25 #define umfpack_UMF_symbolic umfpack_zl_symbolic 26 #define umfpack_UMF_defaults umfpack_zl_defaults 27 28 #else 29 #define umfpack_UMF_free_symbolic umfpack_dl_free_symbolic 30 #define umfpack_UMF_free_numeric umfpack_dl_free_numeric 31 #define umfpack_UMF_wsolve umfpack_dl_wsolve 32 #define umfpack_UMF_numeric umfpack_dl_numeric 33 #define umfpack_UMF_report_numeric umfpack_dl_report_numeric 34 #define umfpack_UMF_report_control umfpack_dl_report_control 35 #define umfpack_UMF_report_status umfpack_dl_report_status 36 #define umfpack_UMF_report_info umfpack_dl_report_info 37 #define umfpack_UMF_report_symbolic umfpack_dl_report_symbolic 38 #define umfpack_UMF_qsymbolic umfpack_dl_qsymbolic 39 #define umfpack_UMF_symbolic umfpack_dl_symbolic 40 #define umfpack_UMF_defaults umfpack_dl_defaults 41 #endif 42 43 #else 44 #if defined(PETSC_USE_COMPLEX) 45 #define umfpack_UMF_free_symbolic umfpack_zi_free_symbolic 46 #define umfpack_UMF_free_numeric umfpack_zi_free_numeric 47 #define umfpack_UMF_wsolve umfpack_zi_wsolve 48 #define umfpack_UMF_numeric umfpack_zi_numeric 49 #define umfpack_UMF_report_numeric umfpack_zi_report_numeric 50 #define umfpack_UMF_report_control umfpack_zi_report_control 51 #define umfpack_UMF_report_status umfpack_zi_report_status 52 #define umfpack_UMF_report_info umfpack_zi_report_info 53 #define umfpack_UMF_report_symbolic umfpack_zi_report_symbolic 54 #define umfpack_UMF_qsymbolic umfpack_zi_qsymbolic 55 #define umfpack_UMF_symbolic umfpack_zi_symbolic 56 #define umfpack_UMF_defaults umfpack_zi_defaults 57 58 #else 59 #define umfpack_UMF_free_symbolic umfpack_di_free_symbolic 60 #define umfpack_UMF_free_numeric umfpack_di_free_numeric 61 #define umfpack_UMF_wsolve umfpack_di_wsolve 62 #define umfpack_UMF_numeric umfpack_di_numeric 63 #define umfpack_UMF_report_numeric umfpack_di_report_numeric 64 #define umfpack_UMF_report_control umfpack_di_report_control 65 #define umfpack_UMF_report_status umfpack_di_report_status 66 #define umfpack_UMF_report_info umfpack_di_report_info 67 #define umfpack_UMF_report_symbolic umfpack_di_report_symbolic 68 #define umfpack_UMF_qsymbolic umfpack_di_qsymbolic 69 #define umfpack_UMF_symbolic umfpack_di_symbolic 70 #define umfpack_UMF_defaults umfpack_di_defaults 71 #endif 72 #endif 73 74 75 #define UF_long long long 76 #define UF_long_max LONG_LONG_MAX 77 #define UF_long_id "%lld" 78 79 EXTERN_C_BEGIN 80 #include <umfpack.h> 81 EXTERN_C_END 82 83 static const char *const UmfpackOrderingTypes[] = {"CHOLMOD","AMD","GIVEN","METIS","BEST","NONE","USER","UmfpackOrderingTypes","UMFPACK_ORDERING_",0}; 84 85 typedef struct { 86 void *Symbolic, *Numeric; 87 double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W; 88 PetscInt *Wi,*perm_c; 89 Mat A; /* Matrix used for factorization */ 90 MatStructure flg; 91 PetscBool PetscMatOrdering; 92 93 /* Flag to clean up UMFPACK objects during Destroy */ 94 PetscBool CleanUpUMFPACK; 95 } Mat_UMFPACK; 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "MatDestroy_UMFPACK" 99 static PetscErrorCode MatDestroy_UMFPACK(Mat A) 100 { 101 PetscErrorCode ierr; 102 Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr; 103 104 PetscFunctionBegin; 105 if (lu && lu->CleanUpUMFPACK) { 106 umfpack_UMF_free_symbolic(&lu->Symbolic); 107 umfpack_UMF_free_numeric(&lu->Numeric); 108 ierr = PetscFree(lu->Wi);CHKERRQ(ierr); 109 ierr = PetscFree(lu->W);CHKERRQ(ierr); 110 ierr = PetscFree(lu->perm_c);CHKERRQ(ierr); 111 } 112 ierr = MatDestroy(&lu->A);CHKERRQ(ierr); 113 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 114 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "MatSolve_UMFPACK_Private" 120 static PetscErrorCode MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag) 121 { 122 Mat_UMFPACK *lu = (Mat_UMFPACK*)A->spptr; 123 Mat_SeqAIJ *a = (Mat_SeqAIJ*)lu->A->data; 124 PetscScalar *av = a->a,*xa; 125 const PetscScalar *ba; 126 PetscErrorCode ierr; 127 PetscInt *ai = a->i,*aj = a->j,status; 128 129 PetscFunctionBegin; 130 /* solve Ax = b by umfpack_*_wsolve */ 131 /* ----------------------------------*/ 132 133 if (!lu->Wi) { /* first time, allocate working space for wsolve */ 134 ierr = PetscMalloc1(A->rmap->n,&lu->Wi);CHKERRQ(ierr); 135 ierr = PetscMalloc1(5*A->rmap->n,&lu->W);CHKERRQ(ierr); 136 } 137 138 ierr = VecGetArrayRead(b,&ba); 139 ierr = VecGetArray(x,&xa); 140 #if defined(PETSC_USE_COMPLEX) 141 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); 142 #else 143 status = umfpack_UMF_wsolve(uflag,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 144 #endif 145 umfpack_UMF_report_info(lu->Control, lu->Info); 146 if (status < 0) { 147 umfpack_UMF_report_status(lu->Control, status); 148 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_wsolve failed"); 149 } 150 151 ierr = VecRestoreArrayRead(b,&ba);CHKERRQ(ierr); 152 ierr = VecRestoreArray(x,&xa);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 #undef __FUNCT__ 157 #define __FUNCT__ "MatSolve_UMFPACK" 158 static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x) 159 { 160 PetscErrorCode ierr; 161 162 PetscFunctionBegin; 163 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 164 ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_Aat);CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 168 #undef __FUNCT__ 169 #define __FUNCT__ "MatSolveTranspose_UMFPACK" 170 static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x) 171 { 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 176 ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "MatLUFactorNumeric_UMFPACK" 182 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info) 183 { 184 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F)->spptr; 185 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 186 PetscInt *ai = a->i,*aj=a->j,status; 187 PetscScalar *av = a->a; 188 PetscErrorCode ierr; 189 190 PetscFunctionBegin; 191 /* numeric factorization of A' */ 192 /* ----------------------------*/ 193 194 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) { 195 umfpack_UMF_free_numeric(&lu->Numeric); 196 } 197 #if defined(PETSC_USE_COMPLEX) 198 status = umfpack_UMF_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 199 #else 200 status = umfpack_UMF_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 201 #endif 202 if (status < 0) { 203 umfpack_UMF_report_status(lu->Control, status); 204 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_numeric failed"); 205 } 206 /* report numeric factorization of A' when Control[PRL] > 3 */ 207 (void) umfpack_UMF_report_numeric(lu->Numeric, lu->Control); 208 209 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 210 ierr = MatDestroy(&lu->A);CHKERRQ(ierr); 211 212 lu->A = A; 213 lu->flg = SAME_NONZERO_PATTERN; 214 lu->CleanUpUMFPACK = PETSC_TRUE; 215 F->ops->solve = MatSolve_UMFPACK; 216 F->ops->solvetranspose = MatSolveTranspose_UMFPACK; 217 PetscFunctionReturn(0); 218 } 219 220 /* 221 Note the r permutation is ignored 222 */ 223 #undef __FUNCT__ 224 #define __FUNCT__ "MatLUFactorSymbolic_UMFPACK" 225 static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 226 { 227 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 228 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F->spptr); 229 PetscErrorCode ierr; 230 PetscInt i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n; 231 PetscScalar *av = a->a; 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 ./configure --download-suitesparse to install PETSc to use UMFPACK 353 354 Consult UMFPACK documentation for more information about the Control parameters 355 which correspond to the options database keys below. 356 357 Options Database Keys: 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 = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 399 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 400 ierr = PetscNewLog(B,&lu);CHKERRQ(ierr); 401 402 B->spptr = lu; 403 B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK; 404 B->ops->destroy = MatDestroy_UMFPACK; 405 B->ops->view = MatView_UMFPACK; 406 407 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_umfpack);CHKERRQ(ierr); 408 409 B->factortype = MAT_FACTOR_LU; 410 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 411 B->preallocated = PETSC_TRUE; 412 413 /* initializations */ 414 /* ------------------------------------------------*/ 415 /* get the default control parameters */ 416 umfpack_UMF_defaults(lu->Control); 417 lu->perm_c = NULL; /* use defaul UMFPACK col permutation */ 418 lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */ 419 420 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr); 421 /* Control parameters used by reporting routiones */ 422 ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],NULL);CHKERRQ(ierr); 423 424 /* Control parameters for symbolic factorization */ 425 ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,3,strategy[0],&idx,&flg);CHKERRQ(ierr); 426 if (flg) { 427 switch (idx) { 428 case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break; 429 case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break; 430 case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break; 431 } 432 } 433 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); 434 if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx; 435 ierr = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],NULL);CHKERRQ(ierr); 436 ierr = PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],NULL);CHKERRQ(ierr); 437 ierr = PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],NULL);CHKERRQ(ierr); 438 ierr = PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],NULL);CHKERRQ(ierr); 439 ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],NULL);CHKERRQ(ierr); 440 ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],NULL);CHKERRQ(ierr); 441 442 /* Control parameters used by numeric factorization */ 443 ierr = PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],NULL);CHKERRQ(ierr); 444 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); 445 ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr); 446 if (flg) { 447 switch (idx) { 448 case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break; 449 case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break; 450 case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break; 451 } 452 } 453 ierr = PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL);CHKERRQ(ierr); 454 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); 455 ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],NULL);CHKERRQ(ierr); 456 457 /* Control parameters used by solve */ 458 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],NULL);CHKERRQ(ierr); 459 460 /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */ 461 ierr = PetscOptionsHasName(NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOrdering);CHKERRQ(ierr); 462 PetscOptionsEnd(); 463 *F = B; 464 PetscFunctionReturn(0); 465 } 466 467 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 468 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 469 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat,MatFactorType,Mat*); 470 471 #undef __FUNCT__ 472 #define __FUNCT__ "MatSolverPackageRegister_SuiteSparse" 473 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_SuiteSparse(void) 474 { 475 PetscErrorCode ierr; 476 477 PetscFunctionBegin; 478 ierr = MatSolverPackageRegister(MATSOLVERUMFPACK,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 479 ierr = MatSolverPackageRegister(MATSOLVERCHOLMOD,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 480 ierr = MatSolverPackageRegister(MATSOLVERCHOLMOD,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 481 ierr = MatSolverPackageRegister(MATSOLVERKLU,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_klu);CHKERRQ(ierr); 482 PetscFunctionReturn(0); 483 } 484