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