1 2 #include "src/ksp/pc/pcimpl.h" /*I "petscpc.h" I*/ 3 #include "src/ksp/pc/impls/factor/factor.h" 4 5 /* Options Database Keys: ??? 6 . -pc_ilu_damping - add damping to diagonal to prevent zero (or very small) pivots 7 . -pc_ilu_shift - apply Manteuffel shift to diagonal to force positive definite preconditioner 8 . -pc_ilu_zeropivot <tol> - set tolerance for what is considered a zero pivot 9 */ 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "PCFactorSetZeroPivot" 13 /*@ 14 PCFactorSetZeroPivot - Sets the size at which smaller pivots are declared to be zero 15 16 Collective on PC 17 18 Input Parameters: 19 + zero - all pivots smaller than this will be considered zero 20 - info - options for factorization 21 22 Options Database Key: 23 . -pc_factor_zeropivot <zero> - Sets tolerance for what is considered a zero pivot 24 25 Level: intermediate 26 27 .keywords: PC, set, factorization, direct, fill 28 29 .seealso: PCFactorSetShiftNonzero(), PCFactorSetShiftPd() 30 @*/ 31 PetscErrorCode PCFactorSetZeroPivot(PetscReal zero,MatFactorInfo *info) 32 { 33 PetscFunctionBegin; 34 info->zeropivot = zero; 35 PetscFunctionReturn(0); 36 } 37 38 #undef __FUNCT__ 39 #define __FUNCT__ "PCFactorSetShiftNonzero" 40 /*@ 41 PCFactorSetShiftNonzero - adds this quantity to the diagonal of the matrix during 42 numerical factorization, thus the matrix has nonzero pivots 43 44 Collective on PC 45 46 Input Parameters: 47 + shift - amount of shift 48 - info - options for factorization 49 50 Options Database Key: 51 . -pc_factor_shiftnonzero <shift> - Sets shift amount or PETSC_DECIDE for the default 52 53 Note: If 0.0 is given, then no shift is used. If a diagonal element is classified as a zero 54 pivot, then the shift is doubled until this is alleviated. 55 56 Level: intermediate 57 58 .keywords: PC, set, factorization, direct, fill 59 60 .seealso: PCFactorSetZeroPivot(), PCFactorSetShiftPd() 61 @*/ 62 PetscErrorCode PCFactorSetShiftNonzero(PetscReal shift,MatFactorInfo *info) 63 { 64 PetscFunctionBegin; 65 if (shift == (PetscReal) PETSC_DECIDE) { 66 info->shiftnz = 1.e-12; 67 } else { 68 info->shiftnz = shift; 69 } 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "PCFactorSetShiftPd" 75 /*@ 76 PCFactorSetShiftPd - specify whether to use Manteuffel shifting. 77 If a matrix factorisation breaks down because of nonpositive pivots, 78 adding sufficient identity to the diagonal will remedy this. 79 Setting this causes a bisection method to find the minimum shift that 80 will lead to a well-defined matrix factor. 81 82 Collective on PC 83 84 Input parameters: 85 + shifting - PETSC_TRUE to set shift else PETSC_FALSE 86 - info - options for factorization 87 88 Options Database Key: 89 . -pc_factor_shiftpd [1/0] - Activate/Deactivate PCFactorSetShiftPd(); the value 90 is optional with 1 being the default 91 92 Level: intermediate 93 94 .keywords: PC, indefinite, factorization 95 96 .seealso: PCFactorSetZeroPivot(), PCFactorSetShiftNonzero() 97 @*/ 98 PetscErrorCode PCFactorSetShiftPd(PetscTruth shifting,MatFactorInfo *info) 99 { 100 PetscFunctionBegin; 101 info->shiftpd = shifting; 102 PetscFunctionReturn(0); 103 } 104 105 /* shift the diagonals when zero pivot is detected */ 106 #undef __FUNCT__ 107 #define __FUNCT__ "PCLUFactorCheckShift" 108 PetscErrorCode PCLUFactorCheckShift(Mat A,MatFactorInfo *info,Mat *B,Shift_Ctx *sctx,PetscInt *newshift) 109 { 110 PetscReal rs; 111 PetscScalar pv; 112 113 PetscFunctionBegin; 114 rs = sctx->rs; 115 pv = sctx->pv; 116 /* printf(" CheckShift: rs: %g\n",rs); */ 117 118 if (PetscAbsScalar(pv) <= info->zeropivot*rs && info->shiftnz){ 119 /* force |diag(*B)| > zeropivot*rs */ 120 if (!sctx->nshift){ 121 sctx->shift_amount = info->shiftnz; 122 } else { 123 sctx->shift_amount *= 2.0; 124 } 125 sctx->lushift = 1; 126 (sctx->nshift)++; 127 *newshift = PETSC_TRUE; 128 } else if (PetscRealPart(pv) <= info->zeropivot*rs && info->shiftpd){ 129 /* force *B to be diagonally dominant */ 130 if (sctx->nshift > sctx->nshift_max) { 131 SETERRQ(PETSC_ERR_CONV_FAILED,"Unable to determine shift to enforce positive definite preconditioner"); 132 } else if (sctx->nshift == sctx->nshift_max) { 133 info->shift_fraction = sctx->shift_hi; 134 sctx->lushift = PETSC_FALSE; 135 } else { 136 sctx->shift_lo = info->shift_fraction; 137 info->shift_fraction = (sctx->shift_hi+sctx->shift_lo)/2.; 138 sctx->lushift = PETSC_TRUE; 139 } 140 sctx->shift_amount = info->shift_fraction * sctx->shift_top; 141 sctx->nshift++; 142 *newshift = PETSC_TRUE; 143 } else if (PetscAbsScalar(pv) <= info->zeropivot*rs){ 144 *newshift = -1; 145 } 146 PetscFunctionReturn(0); 147 } 148 149 150