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