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