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) { 92 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col)); 93 } 94 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 95 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col)); 96 } 97 PetscCall(MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info)); 98 PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info)); 99 dir->hdr.actualfill = info.fill_ratio_needed; 100 } else if (pc->flag != SAME_NONZERO_PATTERN) { 101 PetscBool canuseordering; 102 if (!dir->hdr.reuseordering) { 103 PetscCall(MatDestroy(&((PC_Factor*)dir)->fact)); 104 PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact)); 105 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 106 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering)); 107 if (canuseordering) { 108 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 109 PetscCall(ISDestroy(&dir->col)); 110 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 111 PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col)); 112 if (dir->nonzerosalongdiagonal) { 113 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col)); 114 } 115 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 116 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col)); 117 } 118 } 119 PetscCall(MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info)); 120 PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info)); 121 dir->hdr.actualfill = info.fill_ratio_needed; 122 } else { 123 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 124 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) { 125 PetscCall(MatFactorClearError(((PC_Factor*)dir)->fact)); 126 pc->failedreason = PC_NOERROR; 127 } 128 } 129 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 130 if (err) { /* FactorSymbolic() fails */ 131 pc->failedreason = (PCFailedReason)err; 132 PetscFunctionReturn(0); 133 } 134 135 PetscCall(MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info)); 136 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 137 if (err) { /* FactorNumeric() fails */ 138 pc->failedreason = (PCFailedReason)err; 139 } 140 141 } 142 143 PetscCall(PCFactorGetMatSolverType(pc,&stype)); 144 if (!stype) { 145 MatSolverType solverpackage; 146 PetscCall(MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage)); 147 PetscCall(PCFactorSetMatSolverType(pc,solverpackage)); 148 } 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode PCReset_LU(PC pc) 153 { 154 PC_LU *dir = (PC_LU*)pc->data; 155 156 PetscFunctionBegin; 157 if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) PetscCall(MatDestroy(&((PC_Factor*)dir)->fact)); 158 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row)); 159 PetscCall(ISDestroy(&dir->col)); 160 PetscFunctionReturn(0); 161 } 162 163 static PetscErrorCode PCDestroy_LU(PC pc) 164 { 165 PC_LU *dir = (PC_LU*)pc->data; 166 167 PetscFunctionBegin; 168 PetscCall(PCReset_LU(pc)); 169 PetscCall(PetscFree(((PC_Factor*)dir)->ordering)); 170 PetscCall(PetscFree(((PC_Factor*)dir)->solvertype)); 171 PetscCall(PCFactorClearComposedFunctions(pc)); 172 PetscCall(PetscFree(pc->data)); 173 PetscFunctionReturn(0); 174 } 175 176 static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y) 177 { 178 PC_LU *dir = (PC_LU*)pc->data; 179 180 PetscFunctionBegin; 181 if (dir->hdr.inplace) { 182 PetscCall(MatSolve(pc->pmat,x,y)); 183 } else { 184 PetscCall(MatSolve(((PC_Factor*)dir)->fact,x,y)); 185 } 186 PetscFunctionReturn(0); 187 } 188 189 static PetscErrorCode PCMatApply_LU(PC pc,Mat X,Mat Y) 190 { 191 PC_LU *dir = (PC_LU*)pc->data; 192 193 PetscFunctionBegin; 194 if (dir->hdr.inplace) { 195 PetscCall(MatMatSolve(pc->pmat,X,Y)); 196 } else { 197 PetscCall(MatMatSolve(((PC_Factor*)dir)->fact,X,Y)); 198 } 199 PetscFunctionReturn(0); 200 } 201 202 static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y) 203 { 204 PC_LU *dir = (PC_LU*)pc->data; 205 206 PetscFunctionBegin; 207 if (dir->hdr.inplace) { 208 PetscCall(MatSolveTranspose(pc->pmat,x,y)); 209 } else { 210 PetscCall(MatSolveTranspose(((PC_Factor*)dir)->fact,x,y)); 211 } 212 PetscFunctionReturn(0); 213 } 214 215 /* -----------------------------------------------------------------------------------*/ 216 217 /*MC 218 PCLU - Uses a direct solver, based on LU factorization, as a preconditioner 219 220 Options Database Keys: 221 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 222 . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu 223 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 224 . -pc_factor_fill <fill> - Sets fill amount 225 . -pc_factor_in_place - Activates in-place factorization 226 . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 227 . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase 228 stability of factorization. 229 . -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types 230 . -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default 231 . -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the diagonal. 232 - -mat_solvertype_optionname - options for a specific solver package, for example -mat_mumps_cntl_1 233 234 Notes: 235 Not all options work for all matrix formats 236 Run with -help to see additional options for particular matrix formats or factorization 237 algorithms 238 239 Level: beginner 240 241 Notes: 242 Usually this will compute an "exact" solution in one iteration and does 243 not need a Krylov method (i.e. you can use -ksp_type preonly, or 244 KSPSetType(ksp,KSPPREONLY) for the Krylov method 245 246 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 247 `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 248 `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`, 249 `PCFactorSetPivotInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 250 `PCFactorReorderForNonzeroDiagonal()` 251 M*/ 252 253 PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc) 254 { 255 PC_LU *dir; 256 257 PetscFunctionBegin; 258 PetscCall(PetscNewLog(pc,&dir)); 259 pc->data = (void*)dir; 260 PetscCall(PCFactorInitialize(pc,MAT_FACTOR_LU)); 261 dir->nonzerosalongdiagonal = PETSC_FALSE; 262 263 ((PC_Factor*)dir)->info.fill = 5.0; 264 ((PC_Factor*)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */ 265 ((PC_Factor*)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE; 266 dir->col = NULL; 267 dir->row = NULL; 268 269 pc->ops->reset = PCReset_LU; 270 pc->ops->destroy = PCDestroy_LU; 271 pc->ops->apply = PCApply_LU; 272 pc->ops->matapply = PCMatApply_LU; 273 pc->ops->applytranspose = PCApplyTranspose_LU; 274 pc->ops->setup = PCSetUp_LU; 275 pc->ops->setfromoptions = PCSetFromOptions_LU; 276 pc->ops->view = PCView_Factor; 277 pc->ops->applyrichardson = NULL; 278 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU)); 279 PetscFunctionReturn(0); 280 } 281