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 #include <../src/ksp/pc/impls/factor/factor.h> /*I "petscpc.h" I*/ 8 9 typedef struct { 10 PC_Factor hdr; 11 IS row,col; /* index sets used for reordering */ 12 } PC_Cholesky; 13 14 static PetscErrorCode PCSetFromOptions_Cholesky(PetscOptionItems *PetscOptionsObject,PC pc) 15 { 16 PetscFunctionBegin; 17 PetscOptionsHeadBegin(PetscOptionsObject,"Cholesky options"); 18 PetscCall(PCSetFromOptions_Factor(PetscOptionsObject,pc)); 19 PetscOptionsHeadEnd(); 20 PetscFunctionReturn(0); 21 } 22 23 static PetscErrorCode PCSetUp_Cholesky(PC pc) 24 { 25 PetscBool flg; 26 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 27 MatSolverType stype; 28 MatFactorError err; 29 const char *prefix; 30 31 PetscFunctionBegin; 32 if (!((PetscObject)pc->pmat)->prefix) { 33 PetscCall(PCGetOptionsPrefix(pc,&prefix)); 34 PetscCall(MatSetOptionsPrefix(pc->pmat,prefix)); 35 } 36 pc->failedreason = PC_NOERROR; 37 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill; 38 39 PetscCall(MatSetErrorIfFailure(pc->pmat,pc->erroriffailure)); 40 if (dir->hdr.inplace) { 41 if (dir->row && dir->col && (dir->row != dir->col)) { 42 PetscCall(ISDestroy(&dir->row)); 43 } 44 PetscCall(ISDestroy(&dir->col)); 45 /* should only get reordering if the factor matrix uses it but cannot determine because MatGetFactor() not called */ 46 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 47 PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col)); 48 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 49 PetscCall(ISDestroy(&dir->col)); 50 } 51 if (dir->row) PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 52 PetscCall(MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info)); 53 PetscCall(MatFactorGetError(pc->pmat,&err)); 54 if (err) { /* Factor() fails */ 55 pc->failedreason = (PCFailedReason)err; 56 PetscFunctionReturn(0); 57 } 58 59 ((PC_Factor*)dir)->fact = pc->pmat; 60 } else { 61 MatInfo info; 62 63 if (!pc->setupcalled) { 64 PetscBool canuseordering; 65 if (!((PC_Factor*)dir)->fact) { 66 PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact)); 67 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 68 } 69 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering)); 70 if (canuseordering) { 71 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 72 PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col)); 73 /* check if dir->row == dir->col */ 74 if (dir->row) { 75 PetscCall(ISEqual(dir->row,dir->col,&flg)); 76 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must be equal"); 77 } 78 PetscCall(ISDestroy(&dir->col)); /* only pass one ordering into CholeskyFactor */ 79 80 flg = PETSC_FALSE; 81 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL)); 82 if (flg) { 83 PetscReal tol = 1.e-10; 84 PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL)); 85 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row)); 86 } 87 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 88 } 89 PetscCall(MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info)); 90 PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info)); 91 dir->hdr.actualfill = info.fill_ratio_needed; 92 } else if (pc->flag != SAME_NONZERO_PATTERN) { 93 if (!dir->hdr.reuseordering) { 94 PetscBool canuseordering; 95 PetscCall(MatDestroy(&((PC_Factor*)dir)->fact)); 96 PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact)); 97 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 98 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering)); 99 if (canuseordering) { 100 PetscCall(ISDestroy(&dir->row)); 101 PetscCall(PCFactorSetDefaultOrdering_Factor(pc)); 102 PetscCall(MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col)); 103 PetscCall(ISDestroy(&dir->col)); /* only use dir->row ordering in CholeskyFactor */ 104 105 flg = PETSC_FALSE; 106 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL)); 107 if (flg) { 108 PetscReal tol = 1.e-10; 109 PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL)); 110 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row)); 111 } 112 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row)); 113 } 114 } 115 PetscCall(MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info)); 116 PetscCall(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info)); 117 dir->hdr.actualfill = info.fill_ratio_needed; 118 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 119 } else { 120 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 121 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) { 122 PetscCall(MatFactorClearError(((PC_Factor*)dir)->fact)); 123 pc->failedreason = PC_NOERROR; 124 } 125 } 126 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 127 if (err) { /* FactorSymbolic() fails */ 128 pc->failedreason = (PCFailedReason)err; 129 PetscFunctionReturn(0); 130 } 131 132 PetscCall(MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info)); 133 PetscCall(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 134 if (err) { /* FactorNumeric() fails */ 135 pc->failedreason = (PCFailedReason)err; 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_Cholesky(PC pc) 149 { 150 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 151 152 PetscFunctionBegin; 153 if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) PetscCall(MatDestroy(&((PC_Factor*)dir)->fact)); 154 PetscCall(ISDestroy(&dir->row)); 155 PetscCall(ISDestroy(&dir->col)); 156 PetscFunctionReturn(0); 157 } 158 159 static PetscErrorCode PCDestroy_Cholesky(PC pc) 160 { 161 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 162 163 PetscFunctionBegin; 164 PetscCall(PCReset_Cholesky(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_Cholesky(PC pc,Vec x,Vec y) 173 { 174 PC_Cholesky *dir = (PC_Cholesky*)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_Cholesky(PC pc,Mat X,Mat Y) 186 { 187 PC_Cholesky *dir = (PC_Cholesky*)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 PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y) 199 { 200 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 201 202 PetscFunctionBegin; 203 if (dir->hdr.inplace) { 204 PetscCall(MatForwardSolve(pc->pmat,x,y)); 205 } else { 206 PetscCall(MatForwardSolve(((PC_Factor*)dir)->fact,x,y)); 207 } 208 PetscFunctionReturn(0); 209 } 210 211 static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y) 212 { 213 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 214 215 PetscFunctionBegin; 216 if (dir->hdr.inplace) { 217 PetscCall(MatBackwardSolve(pc->pmat,x,y)); 218 } else { 219 PetscCall(MatBackwardSolve(((PC_Factor*)dir)->fact,x,y)); 220 } 221 PetscFunctionReturn(0); 222 } 223 224 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y) 225 { 226 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 227 228 PetscFunctionBegin; 229 if (dir->hdr.inplace) { 230 PetscCall(MatSolveTranspose(pc->pmat,x,y)); 231 } else { 232 PetscCall(MatSolveTranspose(((PC_Factor*)dir)->fact,x,y)); 233 } 234 PetscFunctionReturn(0); 235 } 236 237 /* -----------------------------------------------------------------------------------*/ 238 239 /* -----------------------------------------------------------------------------------*/ 240 241 /*@ 242 PCFactorSetReuseOrdering - When similar matrices are factored, this 243 causes the ordering computed in the first factor to be used for all 244 following factors. 245 246 Logically Collective on PC 247 248 Input Parameters: 249 + pc - the preconditioner context 250 - flag - PETSC_TRUE to reuse else PETSC_FALSE 251 252 Options Database Key: 253 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 254 255 Level: intermediate 256 257 .seealso: `PCFactorSetReuseFill()` 258 @*/ 259 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 260 { 261 PetscFunctionBegin; 262 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 263 PetscValidLogicalCollectiveBool(pc,flag,2); 264 PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag)); 265 PetscFunctionReturn(0); 266 } 267 268 /*MC 269 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 270 271 Options Database Keys: 272 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 273 . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu 274 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 275 . -pc_factor_fill <fill> - Sets fill amount 276 . -pc_factor_in_place - Activates in-place factorization 277 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 278 279 Notes: 280 Not all options work for all matrix formats 281 282 Level: beginner 283 284 Notes: 285 Usually this will compute an "exact" solution in one iteration and does 286 not need a Krylov method (i.e. you can use -ksp_type preonly, or 287 KSPSetType(ksp,KSPPREONLY) for the Krylov method 288 289 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 290 `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 291 `PCFactorSetFill()`, `PCFactorSetShiftNonzero()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 292 `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()` 293 294 M*/ 295 296 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) 297 { 298 PC_Cholesky *dir; 299 300 PetscFunctionBegin; 301 PetscCall(PetscNewLog(pc,&dir)); 302 pc->data = (void*)dir; 303 PetscCall(PCFactorInitialize(pc,MAT_FACTOR_CHOLESKY)); 304 305 ((PC_Factor*)dir)->info.fill = 5.0; 306 307 pc->ops->destroy = PCDestroy_Cholesky; 308 pc->ops->reset = PCReset_Cholesky; 309 pc->ops->apply = PCApply_Cholesky; 310 pc->ops->matapply = PCMatApply_Cholesky; 311 pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky; 312 pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky; 313 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 314 pc->ops->setup = PCSetUp_Cholesky; 315 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 316 pc->ops->view = PCView_Factor; 317 pc->ops->applyrichardson = NULL; 318 PetscFunctionReturn(0); 319 } 320