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 PetscReal actualfill; /* actual fill in factor */ 12 PetscBool inplace; /* flag indicating in-place factorization */ 13 IS row,col; /* index sets used for reordering */ 14 PetscBool reuseordering; /* reuses previous reordering computed */ 15 PetscBool reusefill; /* reuse fill from previous Cholesky */ 16 } PC_Cholesky; 17 18 #undef __FUNCT__ 19 #define __FUNCT__ "PCFactorSetReuseOrdering_Cholesky" 20 static PetscErrorCode PCFactorSetReuseOrdering_Cholesky(PC pc,PetscBool flag) 21 { 22 PC_Cholesky *lu = (PC_Cholesky*)pc->data; 23 24 PetscFunctionBegin; 25 lu->reuseordering = flag; 26 PetscFunctionReturn(0); 27 } 28 29 #undef __FUNCT__ 30 #define __FUNCT__ "PCFactorSetReuseFill_Cholesky" 31 static PetscErrorCode PCFactorSetReuseFill_Cholesky(PC pc,PetscBool flag) 32 { 33 PC_Cholesky *lu = (PC_Cholesky*)pc->data; 34 35 PetscFunctionBegin; 36 lu->reusefill = flag; 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "PCSetFromOptions_Cholesky" 42 static PetscErrorCode PCSetFromOptions_Cholesky(PetscOptions *PetscOptionsObject,PC pc) 43 { 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = PetscOptionsHead(PetscOptionsObject,"Cholesky options");CHKERRQ(ierr); 48 ierr = PCSetFromOptions_Factor(PetscOptionsObject,pc);CHKERRQ(ierr); 49 ierr = PetscOptionsTail();CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 #undef __FUNCT__ 54 #define __FUNCT__ "PCView_Cholesky" 55 static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer) 56 { 57 PC_Cholesky *chol = (PC_Cholesky*)pc->data; 58 PetscErrorCode ierr; 59 PetscBool iascii; 60 61 PetscFunctionBegin; 62 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 63 if (iascii) { 64 if (chol->inplace) { 65 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");CHKERRQ(ierr); 66 } else { 67 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");CHKERRQ(ierr); 68 } 69 70 if (chol->reusefill) {ierr = PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");CHKERRQ(ierr);} 71 if (chol->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");CHKERRQ(ierr);} 72 } 73 ierr = PCView_Factor(pc,viewer);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "PCSetUp_Cholesky" 80 static PetscErrorCode PCSetUp_Cholesky(PC pc) 81 { 82 PetscErrorCode ierr; 83 PetscBool flg; 84 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 85 86 PetscFunctionBegin; 87 if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill; 88 89 if (dir->inplace) { 90 if (dir->row && dir->col && (dir->row != dir->col)) { 91 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 92 } 93 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 94 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 95 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 96 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 97 } 98 if (dir->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr);} 99 ierr = MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 100 101 ((PC_Factor*)dir)->fact = pc->pmat; 102 } else { 103 MatInfo info; 104 if (!pc->setupcalled) { 105 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 106 /* check if dir->row == dir->col */ 107 ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr); 108 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal"); 109 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */ 110 111 flg = PETSC_FALSE; 112 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);CHKERRQ(ierr); 113 if (flg) { 114 PetscReal tol = 1.e-10; 115 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);CHKERRQ(ierr); 116 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 117 } 118 if (dir->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr);} 119 if (!((PC_Factor*)dir)->fact) { 120 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 121 } 122 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 123 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 124 dir->actualfill = info.fill_ratio_needed; 125 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);CHKERRQ(ierr); 126 } else if (pc->flag != SAME_NONZERO_PATTERN) { 127 if (!dir->reuseordering) { 128 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 129 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 130 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only use dir->row ordering in CholeskyFactor */ 131 132 flg = PETSC_FALSE; 133 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);CHKERRQ(ierr); 134 if (flg) { 135 PetscReal tol = 1.e-10; 136 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);CHKERRQ(ierr); 137 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 138 } 139 if (dir->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr);} 140 } 141 ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 142 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 143 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 144 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 145 dir->actualfill = info.fill_ratio_needed; 146 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);CHKERRQ(ierr); 147 } 148 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "PCReset_Cholesky" 155 static PetscErrorCode PCReset_Cholesky(PC pc) 156 { 157 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 if (!dir->inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 162 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 163 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 164 PetscFunctionReturn(0); 165 } 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "PCDestroy_Cholesky" 169 static PetscErrorCode PCDestroy_Cholesky(PC pc) 170 { 171 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 ierr = PCReset_Cholesky(pc);CHKERRQ(ierr); 176 ierr = PetscFree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 177 ierr = PetscFree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 178 ierr = PetscFree(pc->data);CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "PCApply_Cholesky" 184 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 185 { 186 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 187 PetscErrorCode ierr; 188 189 PetscFunctionBegin; 190 if (dir->inplace) { 191 ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr); 192 } else { 193 ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "PCApplyTranspose_Cholesky" 200 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y) 201 { 202 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 if (dir->inplace) { 207 ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr); 208 } else { 209 ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 210 } 211 PetscFunctionReturn(0); 212 } 213 214 /* -----------------------------------------------------------------------------------*/ 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky" 218 static PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc,PetscBool flg) 219 { 220 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 221 222 PetscFunctionBegin; 223 dir->inplace = flg; 224 PetscFunctionReturn(0); 225 } 226 227 #undef __FUNCT__ 228 #define __FUNCT__ "PCFactorGetUseInPlace_Cholesky" 229 static PetscErrorCode PCFactorGetUseInPlace_Cholesky(PC pc,PetscBool *flg) 230 { 231 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 232 233 PetscFunctionBegin; 234 *flg = dir->inplace; 235 PetscFunctionReturn(0); 236 } 237 238 /* -----------------------------------------------------------------------------------*/ 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "PCFactorSetReuseOrdering" 242 /*@ 243 PCFactorSetReuseOrdering - When similar matrices are factored, this 244 causes the ordering computed in the first factor to be used for all 245 following factors. 246 247 Logically Collective on PC 248 249 Input Parameters: 250 + pc - the preconditioner context 251 - flag - PETSC_TRUE to reuse else PETSC_FALSE 252 253 Options Database Key: 254 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 255 256 Level: intermediate 257 258 .keywords: PC, levels, reordering, factorization, incomplete, LU 259 260 .seealso: PCFactorSetReuseFill() 261 @*/ 262 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 263 { 264 PetscErrorCode ierr; 265 266 PetscFunctionBegin; 267 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 268 PetscValidLogicalCollectiveBool(pc,flag,2); 269 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 273 /*MC 274 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 275 276 Options Database Keys: 277 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 278 . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu 279 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 280 . -pc_factor_fill <fill> - Sets fill amount 281 . -pc_factor_in_place - Activates in-place factorization 282 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 283 284 Notes: Not all options work for all matrix formats 285 286 Level: beginner 287 288 Concepts: Cholesky factorization, direct solver 289 290 Notes: Usually this will compute an "exact" solution in one iteration and does 291 not need a Krylov method (i.e. you can use -ksp_type preonly, or 292 KSPSetType(ksp,KSPPREONLY) for the Krylov method 293 294 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 295 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 296 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 297 PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType() 298 299 M*/ 300 301 #undef __FUNCT__ 302 #define __FUNCT__ "PCCreate_Cholesky" 303 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) 304 { 305 PetscErrorCode ierr; 306 PC_Cholesky *dir; 307 308 PetscFunctionBegin; 309 ierr = PetscNewLog(pc,&dir);CHKERRQ(ierr); 310 311 ((PC_Factor*)dir)->fact = 0; 312 dir->inplace = PETSC_FALSE; 313 314 ierr = MatFactorInfoInitialize(&((PC_Factor*)dir)->info);CHKERRQ(ierr); 315 316 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 317 ((PC_Factor*)dir)->info.fill = 5.0; 318 ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE; 319 ((PC_Factor*)dir)->info.shiftamount = 0.0; 320 ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON; 321 ((PC_Factor*)dir)->info.pivotinblocks = 1.0; 322 323 dir->col = 0; 324 dir->row = 0; 325 326 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 327 ierr = PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 328 329 dir->reusefill = PETSC_FALSE; 330 dir->reuseordering = PETSC_FALSE; 331 pc->data = (void*)dir; 332 333 pc->ops->destroy = PCDestroy_Cholesky; 334 pc->ops->reset = PCReset_Cholesky; 335 pc->ops->apply = PCApply_Cholesky; 336 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 337 pc->ops->setup = PCSetUp_Cholesky; 338 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 339 pc->ops->view = PCView_Cholesky; 340 pc->ops->applyrichardson = 0; 341 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 342 343 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr); 344 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 345 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 346 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 347 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);CHKERRQ(ierr); 348 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 349 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);CHKERRQ(ierr); 350 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr); 351 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_Cholesky);CHKERRQ(ierr); 352 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 353 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr); 354 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr); 355 PetscFunctionReturn(0); 356 } 357