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(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