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