1 2 /* 3 Defines a direct factorization preconditioner for any Mat implementation 4 Note: this need not be consided a preconditioner since it supplies 5 a direct solver. 6 */ 7 8 #include <../src/ksp/pc/impls/factor/lu/lu.h> /*I "petscpc.h" I*/ 9 10 PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc,PetscReal z) 11 { 12 PC_LU *lu = (PC_LU*)pc->data; 13 14 PetscFunctionBegin; 15 lu->nonzerosalongdiagonal = PETSC_TRUE; 16 if (z == PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10; 17 else lu->nonzerosalongdiagonaltol = z; 18 PetscFunctionReturn(0); 19 } 20 21 static PetscErrorCode PCSetFromOptions_LU(PetscOptionItems *PetscOptionsObject,PC pc) 22 { 23 PC_LU *lu = (PC_LU*)pc->data; 24 PetscBool flg = PETSC_FALSE; 25 PetscReal tol; 26 27 PetscFunctionBegin; 28 PetscOptionsHeadBegin(PetscOptionsObject,"LU options"); 29 PetscCall(PCSetFromOptions_Factor(PetscOptionsObject,pc)); 30 31 PetscCall(PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg)); 32 if (flg) { 33 tol = PETSC_DECIDE; 34 PetscCall(PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,NULL)); 35 PetscCall(PCFactorReorderForNonzeroDiagonal(pc,tol)); 36 } 37 PetscOptionsHeadEnd(); 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode PCSetUp_LU(PC pc) 42 { 43 PC_LU *dir = (PC_LU*)pc->data; 44 MatSolverType stype; 45 MatFactorError err; 46 const char *prefix; 47 48 PetscFunctionBegin; 49 pc->failedreason = PC_NOERROR; 50 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill; 51 52 PetscCall(PCGetOptionsPrefix(pc,&prefix)); 53 PetscCall(MatSetOptionsPrefixFactor(pc->pmat,prefix)); 54 55 PetscCall(MatSetErrorIfFailure(pc->pmat,pc->erroriffailure)); 56 if (dir->hdr.inplace) { 57 MatFactorType ftype; 58 59 PetscCall(MatGetFactorType(pc->pmat, &ftype)); 60 if (ftype == MAT_FACTOR_NONE) { 61 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 62 PetscCall(ISDestroy(&dir->col)); 63 /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */ 64 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 65 PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col)); 66 if (dir->row) { 67 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 68 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col)); 69 } 70 PetscCall(MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info)); 71 PetscCall(MatFactorGetError(pc->pmat,&err)); 72 if (err) { /* Factor() fails */ 73 pc->failedreason = (PCFailedReason)err; 74 PetscFunctionReturn(0); 75 } 76 } 77 ((PC_Factor*)dir)->fact = pc->pmat; 78 } else { 79 MatInfo info; 80 81 if (!pc->setupcalled) { 82 PetscBool canuseordering; 83 if (!((PC_Factor*)dir)->fact) { 84 PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact)); 85 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 86 } 87 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering)); 88 if (canuseordering) { 89 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 90 PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col)); 91 if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col)); 92 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 93 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col)); 94 } 95 PetscCall(MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info)); 96 PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info)); 97 dir->hdr.actualfill = info.fill_ratio_needed; 98 } else if (pc->flag != SAME_NONZERO_PATTERN) { 99 PetscBool canuseordering; 100 if (!dir->hdr.reuseordering) { 101 PetscCall(MatDestroy(&((PC_Factor*)dir)->fact)); 102 PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact)); 103 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 104 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering)); 105 if (canuseordering) { 106 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 107 PetscCall(ISDestroy(&dir->col)); 108 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 109 PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col)); 110 if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col)); 111 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 112 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col)); 113 } 114 } 115 PetscCall(MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info)); 116 PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info)); 117 dir->hdr.actualfill = info.fill_ratio_needed; 118 } else { 119 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 120 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) { 121 PetscCall(MatFactorClearError(((PC_Factor*)dir)->fact)); 122 pc->failedreason = PC_NOERROR; 123 } 124 } 125 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 126 if (err) { /* FactorSymbolic() fails */ 127 pc->failedreason = (PCFailedReason)err; 128 PetscFunctionReturn(0); 129 } 130 131 PetscCall(MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info)); 132 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 133 if (err) { /* FactorNumeric() fails */ 134 pc->failedreason = (PCFailedReason)err; 135 } 136 137 } 138 139 PetscCall(PCFactorGetMatSolverType(pc,&stype)); 140 if (!stype) { 141 MatSolverType solverpackage; 142 PetscCall(MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage)); 143 PetscCall(PCFactorSetMatSolverType(pc,solverpackage)); 144 } 145 PetscFunctionReturn(0); 146 } 147 148 static PetscErrorCode PCReset_LU(PC pc) 149 { 150 PC_LU *dir = (PC_LU*)pc->data; 151 152 PetscFunctionBegin; 153 if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) PetscCall(MatDestroy(&((PC_Factor*)dir)->fact)); 154 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 155 PetscCall(ISDestroy(&dir->col)); 156 PetscFunctionReturn(0); 157 } 158 159 static PetscErrorCode PCDestroy_LU(PC pc) 160 { 161 PC_LU *dir = (PC_LU*)pc->data; 162 163 PetscFunctionBegin; 164 PetscCall(PCReset_LU(pc)); 165 PetscCall(PetscFree(((PC_Factor*)dir)->ordering)); 166 PetscCall(PetscFree(((PC_Factor*)dir)->solvertype)); 167 PetscCall(PCFactorClearComposedFunctions(pc)); 168 PetscCall(PetscFree(pc->data)); 169 PetscFunctionReturn(0); 170 } 171 172 static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y) 173 { 174 PC_LU *dir = (PC_LU*)pc->data; 175 176 PetscFunctionBegin; 177 if (dir->hdr.inplace) { 178 PetscCall(MatSolve(pc->pmat,x,y)); 179 } else { 180 PetscCall(MatSolve(((PC_Factor*)dir)->fact,x,y)); 181 } 182 PetscFunctionReturn(0); 183 } 184 185 static PetscErrorCode PCMatApply_LU(PC pc,Mat X,Mat Y) 186 { 187 PC_LU *dir = (PC_LU*)pc->data; 188 189 PetscFunctionBegin; 190 if (dir->hdr.inplace) { 191 PetscCall(MatMatSolve(pc->pmat,X,Y)); 192 } else { 193 PetscCall(MatMatSolve(((PC_Factor*)dir)->fact,X,Y)); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y) 199 { 200 PC_LU *dir = (PC_LU*)pc->data; 201 202 PetscFunctionBegin; 203 if (dir->hdr.inplace) { 204 PetscCall(MatSolveTranspose(pc->pmat,x,y)); 205 } else { 206 PetscCall(MatSolveTranspose(((PC_Factor*)dir)->fact,x,y)); 207 } 208 PetscFunctionReturn(0); 209 } 210 211 /* -----------------------------------------------------------------------------------*/ 212 213 /*MC 214 PCLU - Uses a direct solver, based on LU factorization, as a preconditioner 215 216 Options Database Keys: 217 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 218 . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu 219 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 220 . -pc_factor_fill <fill> - Sets fill amount 221 . -pc_factor_in_place - Activates in-place factorization 222 . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 223 . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase 224 stability of factorization. 225 . -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types 226 . -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default 227 . -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the diagonal. 228 - -mat_solvertype_optionname - options for a specific solver package, for example -mat_mumps_cntl_1 229 230 Notes: 231 Not all options work for all matrix formats 232 Run with -help to see additional options for particular matrix formats or factorization 233 algorithms 234 235 Level: beginner 236 237 Notes: 238 Usually this will compute an "exact" solution in one iteration and does 239 not need a Krylov method (i.e. you can use -ksp_type preonly, or 240 KSPSetType(ksp,KSPPREONLY) for the Krylov method 241 242 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 243 `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 244 `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`, 245 `PCFactorSetPivotInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 246 `PCFactorReorderForNonzeroDiagonal()` 247 M*/ 248 249 PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc) 250 { 251 PC_LU *dir; 252 253 PetscFunctionBegin; 254 PetscCall(PetscNewLog(pc,&dir)); 255 pc->data = (void*)dir; 256 PetscCall(PCFactorInitialize(pc,MAT_FACTOR_LU)); 257 dir->nonzerosalongdiagonal = PETSC_FALSE; 258 259 ((PC_Factor*)dir)->info.fill = 5.0; 260 ((PC_Factor*)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */ 261 ((PC_Factor*)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE; 262 dir->col = NULL; 263 dir->row = NULL; 264 265 pc->ops->reset = PCReset_LU; 266 pc->ops->destroy = PCDestroy_LU; 267 pc->ops->apply = PCApply_LU; 268 pc->ops->matapply = PCMatApply_LU; 269 pc->ops->applytranspose = PCApplyTranspose_LU; 270 pc->ops->setup = PCSetUp_LU; 271 pc->ops->setfromoptions = PCSetFromOptions_LU; 272 pc->ops->view = PCView_Factor; 273 pc->ops->applyrichardson = NULL; 274 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU)); 275 PetscFunctionReturn(0); 276 } 277