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