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