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