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