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 = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 45 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 46 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 47 } 48 if (dir->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr);} 49 ierr = MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 50 ierr = MatFactorGetError(pc->pmat,&err);CHKERRQ(ierr); 51 if (err) { /* Factor() fails */ 52 pc->failedreason = (PCFailedReason)err; 53 PetscFunctionReturn(0); 54 } 55 56 ((PC_Factor*)dir)->fact = pc->pmat; 57 } else { 58 MatInfo info; 59 60 if (!pc->setupcalled) { 61 PetscBool useordering; 62 if (!((PC_Factor*)dir)->fact) { 63 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 64 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);CHKERRQ(ierr); 65 } 66 ierr = MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);CHKERRQ(ierr); 67 if (useordering) { 68 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 69 /* check if dir->row == dir->col */ 70 ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr); 71 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal"); 72 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */ 73 74 flg = PETSC_FALSE; 75 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);CHKERRQ(ierr); 76 if (flg) { 77 PetscReal tol = 1.e-10; 78 ierr = PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);CHKERRQ(ierr); 79 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 80 } 81 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr); 82 } 83 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 84 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 85 dir->hdr.actualfill = info.fill_ratio_needed; 86 } else if (pc->flag != SAME_NONZERO_PATTERN) { 87 if (!dir->hdr.reuseordering) { 88 PetscBool useordering; 89 ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 90 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 91 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);CHKERRQ(ierr); 92 ierr = MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);CHKERRQ(ierr); 93 if (useordering) { 94 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 95 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 96 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only use dir->row ordering in CholeskyFactor */ 97 98 flg = PETSC_FALSE; 99 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);CHKERRQ(ierr); 100 if (flg) { 101 PetscReal tol = 1.e-10; 102 ierr = PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);CHKERRQ(ierr); 103 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 104 } 105 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr); 106 } 107 } 108 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 109 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 110 dir->hdr.actualfill = info.fill_ratio_needed; 111 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);CHKERRQ(ierr); 112 } else { 113 ierr = MatFactorGetError(((PC_Factor*)dir)->fact,&err);CHKERRQ(ierr); 114 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) { 115 ierr = MatFactorClearError(((PC_Factor*)dir)->fact);CHKERRQ(ierr); 116 pc->failedreason = PC_NOERROR; 117 } 118 } 119 ierr = MatFactorGetError(((PC_Factor*)dir)->fact,&err);CHKERRQ(ierr); 120 if (err) { /* FactorSymbolic() fails */ 121 pc->failedreason = (PCFailedReason)err; 122 PetscFunctionReturn(0); 123 } 124 125 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 126 ierr = MatFactorGetError(((PC_Factor*)dir)->fact,&err);CHKERRQ(ierr); 127 if (err) { /* FactorNumeric() fails */ 128 pc->failedreason = (PCFailedReason)err; 129 } 130 } 131 132 ierr = PCFactorGetMatSolverType(pc,&stype);CHKERRQ(ierr); 133 if (!stype) { 134 MatSolverType solverpackage; 135 ierr = MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);CHKERRQ(ierr); 136 ierr = PCFactorSetMatSolverType(pc,solverpackage);CHKERRQ(ierr); 137 } 138 PetscFunctionReturn(0); 139 } 140 141 static PetscErrorCode PCReset_Cholesky(PC pc) 142 { 143 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 148 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 149 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 static PetscErrorCode PCDestroy_Cholesky(PC pc) 154 { 155 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PCReset_Cholesky(pc);CHKERRQ(ierr); 160 ierr = PetscFree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 161 ierr = PetscFree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 162 ierr = PetscFree(pc->data);CHKERRQ(ierr); 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 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 if (dir->hdr.inplace) { 173 ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr); 174 } else { 175 ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 176 } 177 PetscFunctionReturn(0); 178 } 179 180 static PetscErrorCode PCMatApply_Cholesky(PC pc,Mat X,Mat Y) 181 { 182 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 if (dir->hdr.inplace) { 187 ierr = MatMatSolve(pc->pmat,X,Y);CHKERRQ(ierr); 188 } else { 189 ierr = MatMatSolve(((PC_Factor*)dir)->fact,X,Y);CHKERRQ(ierr); 190 } 191 PetscFunctionReturn(0); 192 } 193 194 static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y) 195 { 196 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 if (dir->hdr.inplace) { 201 ierr = MatForwardSolve(pc->pmat,x,y);CHKERRQ(ierr); 202 } else { 203 ierr = MatForwardSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 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 PetscErrorCode ierr; 212 213 PetscFunctionBegin; 214 if (dir->hdr.inplace) { 215 ierr = MatBackwardSolve(pc->pmat,x,y);CHKERRQ(ierr); 216 } else { 217 ierr = MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 218 } 219 PetscFunctionReturn(0); 220 } 221 222 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y) 223 { 224 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 225 PetscErrorCode ierr; 226 227 PetscFunctionBegin; 228 if (dir->hdr.inplace) { 229 ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr); 230 } else { 231 ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 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 PetscErrorCode ierr; 261 262 PetscFunctionBegin; 263 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 264 PetscValidLogicalCollectiveBool(pc,flag,2); 265 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 266 PetscFunctionReturn(0); 267 } 268 269 /*MC 270 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 271 272 Options Database Keys: 273 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 274 . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu 275 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 276 . -pc_factor_fill <fill> - Sets fill amount 277 . -pc_factor_in_place - Activates in-place factorization 278 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 279 280 Notes: 281 Not all options work for all matrix formats 282 283 Level: beginner 284 285 Notes: 286 Usually this will compute an "exact" solution in one iteration and does 287 not need a Krylov method (i.e. you can use -ksp_type preonly, or 288 KSPSetType(ksp,KSPPREONLY) for the Krylov method 289 290 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 291 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 292 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 293 PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType() 294 295 M*/ 296 297 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) 298 { 299 PetscErrorCode ierr; 300 PC_Cholesky *dir; 301 302 PetscFunctionBegin; 303 ierr = PetscNewLog(pc,&dir);CHKERRQ(ierr); 304 pc->data = (void*)dir; 305 ierr = PCFactorInitialize(pc);CHKERRQ(ierr); 306 307 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 308 ((PC_Factor*)dir)->info.fill = 5.0; 309 310 dir->col = NULL; 311 dir->row = NULL; 312 313 /* MATORDERINGNATURAL_OR_ND allows selecting type based on matrix type sbaij or aij */ 314 ierr = PetscStrallocpy(MATORDERINGNATURAL_OR_ND,(char**)&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 315 316 pc->ops->destroy = PCDestroy_Cholesky; 317 pc->ops->reset = PCReset_Cholesky; 318 pc->ops->apply = PCApply_Cholesky; 319 pc->ops->matapply = PCMatApply_Cholesky; 320 pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky; 321 pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky; 322 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 323 pc->ops->setup = PCSetUp_Cholesky; 324 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 325 pc->ops->view = PCView_Factor; 326 pc->ops->applyrichardson = NULL; 327 PetscFunctionReturn(0); 328 } 329