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(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 232 diagonal. 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