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 = PetscObjectTypeCompare((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 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 100 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 101 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 102 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 103 } 104 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 105 ierr = MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 106 ((PC_Factor*)dir)->fact = pc->pmat; 107 } else { 108 MatInfo info; 109 if (!pc->setupcalled) { 110 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 111 /* check if dir->row == dir->col */ 112 ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr); 113 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal"); 114 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */ 115 116 flg = PETSC_FALSE; 117 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 118 if (flg) { 119 PetscReal tol = 1.e-10; 120 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 121 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 122 } 123 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 124 if (!((PC_Factor*)dir)->fact){ 125 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 126 } 127 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 128 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 129 dir->actualfill = info.fill_ratio_needed; 130 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 131 } else if (pc->flag != SAME_NONZERO_PATTERN) { 132 if (!dir->reuseordering) { 133 if (dir->row && dir->col && (dir->row != dir->col)) { 134 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 135 } 136 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 137 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 138 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 139 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 140 } 141 flg = PETSC_FALSE; 142 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 143 if (flg) { 144 PetscReal tol = 1.e-10; 145 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 146 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 147 } 148 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 149 } 150 ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 151 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 152 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 153 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 154 dir->actualfill = info.fill_ratio_needed; 155 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 156 } 157 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "PCReset_Cholesky" 164 static PetscErrorCode PCReset_Cholesky(PC pc) 165 { 166 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 167 PetscErrorCode ierr; 168 169 PetscFunctionBegin; 170 if (!dir->inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 171 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 172 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "PCDestroy_Cholesky" 178 static PetscErrorCode PCDestroy_Cholesky(PC pc) 179 { 180 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 ierr = PCReset_Cholesky(pc);CHKERRQ(ierr); 185 ierr = PetscFree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 186 ierr = PetscFree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 187 ierr = PetscFree(pc->data);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "PCApply_Cholesky" 193 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 194 { 195 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 if (dir->inplace) {ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr);} 200 else {ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "PCApplyTranspose_Cholesky" 206 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y) 207 { 208 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 if (dir->inplace) {ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr);} 213 else {ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 214 PetscFunctionReturn(0); 215 } 216 217 /* -----------------------------------------------------------------------------------*/ 218 219 EXTERN_C_BEGIN 220 #undef __FUNCT__ 221 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky" 222 PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc) 223 { 224 PC_Cholesky *dir; 225 226 PetscFunctionBegin; 227 dir = (PC_Cholesky*)pc->data; 228 dir->inplace = PETSC_TRUE; 229 PetscFunctionReturn(0); 230 } 231 EXTERN_C_END 232 233 /* -----------------------------------------------------------------------------------*/ 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "PCFactorSetReuseOrdering" 237 /*@ 238 PCFactorSetReuseOrdering - When similar matrices are factored, this 239 causes the ordering computed in the first factor to be used for all 240 following factors. 241 242 Logically Collective on PC 243 244 Input Parameters: 245 + pc - the preconditioner context 246 - flag - PETSC_TRUE to reuse else PETSC_FALSE 247 248 Options Database Key: 249 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 250 251 Level: intermediate 252 253 .keywords: PC, levels, reordering, factorization, incomplete, LU 254 255 .seealso: PCFactorSetReuseFill() 256 @*/ 257 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 258 { 259 PetscErrorCode ierr; 260 261 PetscFunctionBegin; 262 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 263 PetscValidLogicalCollectiveBool(pc,flag,2); 264 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 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_package - Actives PCFactorSetMatSolverPackage() 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: Not all options work for all matrix formats 281 282 Level: beginner 283 284 Concepts: Cholesky factorization, direct solver 285 286 Notes: 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(), PCFactorSetMatOrderingType() 294 295 M*/ 296 297 EXTERN_C_BEGIN 298 #undef __FUNCT__ 299 #define __FUNCT__ "PCCreate_Cholesky" 300 PetscErrorCode PCCreate_Cholesky(PC pc) 301 { 302 PetscErrorCode ierr; 303 PC_Cholesky *dir; 304 305 PetscFunctionBegin; 306 ierr = PetscNewLog(pc,PC_Cholesky,&dir);CHKERRQ(ierr); 307 308 ((PC_Factor*)dir)->fact = 0; 309 dir->inplace = PETSC_FALSE; 310 ierr = MatFactorInfoInitialize(&((PC_Factor*)dir)->info);CHKERRQ(ierr); 311 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 312 ((PC_Factor*)dir)->info.fill = 5.0; 313 ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE; 314 ((PC_Factor*)dir)->info.shiftamount = 0.0; 315 ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON; 316 ((PC_Factor*)dir)->info.pivotinblocks = 1.0; 317 dir->col = 0; 318 dir->row = 0; 319 ierr = PetscStrallocpy(MATORDERINGNATURAL,&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 320 ierr = PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 321 dir->reusefill = PETSC_FALSE; 322 dir->reuseordering = PETSC_FALSE; 323 pc->data = (void*)dir; 324 325 pc->ops->destroy = PCDestroy_Cholesky; 326 pc->ops->reset = PCReset_Cholesky; 327 pc->ops->apply = PCApply_Cholesky; 328 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 329 pc->ops->setup = PCSetUp_Cholesky; 330 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 331 pc->ops->view = PCView_Cholesky; 332 pc->ops->applyrichardson = 0; 333 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 334 335 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor", 336 PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr); 337 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 338 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 339 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 340 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 342 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 343 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor", 344 PCFactorSetShiftType_Factor);CHKERRQ(ierr); 345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor", 346 PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 347 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 348 PCFactorSetFill_Factor);CHKERRQ(ierr); 349 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky", 350 PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr); 351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 352 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 353 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky", 354 PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr); 355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky", 356 PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr); 357 PetscFunctionReturn(0); 358 } 359 EXTERN_C_END 360