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