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 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 134 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 135 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only use dir->row ordering in CholeskyFactor */ 136 137 flg = PETSC_FALSE; 138 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 139 if (flg) { 140 PetscReal tol = 1.e-10; 141 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 142 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 143 } 144 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 145 } 146 ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 147 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 148 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 149 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 150 dir->actualfill = info.fill_ratio_needed; 151 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 152 } 153 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 154 } 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "PCReset_Cholesky" 160 static PetscErrorCode PCReset_Cholesky(PC pc) 161 { 162 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 163 PetscErrorCode ierr; 164 165 PetscFunctionBegin; 166 if (!dir->inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 167 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 168 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCDestroy_Cholesky" 174 static PetscErrorCode PCDestroy_Cholesky(PC pc) 175 { 176 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 177 PetscErrorCode ierr; 178 179 PetscFunctionBegin; 180 ierr = PCReset_Cholesky(pc);CHKERRQ(ierr); 181 ierr = PetscFree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 182 ierr = PetscFree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 183 ierr = PetscFree(pc->data);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "PCApply_Cholesky" 189 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 190 { 191 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 if (dir->inplace) {ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr);} 196 else {ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 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) {ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr);} 209 else {ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 210 PetscFunctionReturn(0); 211 } 212 213 /* -----------------------------------------------------------------------------------*/ 214 215 EXTERN_C_BEGIN 216 #undef __FUNCT__ 217 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky" 218 PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc) 219 { 220 PC_Cholesky *dir; 221 222 PetscFunctionBegin; 223 dir = (PC_Cholesky*)pc->data; 224 dir->inplace = PETSC_TRUE; 225 PetscFunctionReturn(0); 226 } 227 EXTERN_C_END 228 229 /* -----------------------------------------------------------------------------------*/ 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "PCFactorSetReuseOrdering" 233 /*@ 234 PCFactorSetReuseOrdering - When similar matrices are factored, this 235 causes the ordering computed in the first factor to be used for all 236 following factors. 237 238 Logically Collective on PC 239 240 Input Parameters: 241 + pc - the preconditioner context 242 - flag - PETSC_TRUE to reuse else PETSC_FALSE 243 244 Options Database Key: 245 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 246 247 Level: intermediate 248 249 .keywords: PC, levels, reordering, factorization, incomplete, LU 250 251 .seealso: PCFactorSetReuseFill() 252 @*/ 253 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 254 { 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 259 PetscValidLogicalCollectiveBool(pc,flag,2); 260 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 265 /*MC 266 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 267 268 Options Database Keys: 269 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 270 . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu 271 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 272 . -pc_factor_fill <fill> - Sets fill amount 273 . -pc_factor_in_place - Activates in-place factorization 274 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 275 276 Notes: Not all options work for all matrix formats 277 278 Level: beginner 279 280 Concepts: Cholesky factorization, direct solver 281 282 Notes: Usually this will compute an "exact" solution in one iteration and does 283 not need a Krylov method (i.e. you can use -ksp_type preonly, or 284 KSPSetType(ksp,KSPPREONLY) for the Krylov method 285 286 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 287 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 288 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 289 PCFactorSetUseInPlace(), PCFactorSetMatOrderingType() 290 291 M*/ 292 293 EXTERN_C_BEGIN 294 #undef __FUNCT__ 295 #define __FUNCT__ "PCCreate_Cholesky" 296 PetscErrorCode PCCreate_Cholesky(PC pc) 297 { 298 PetscErrorCode ierr; 299 PC_Cholesky *dir; 300 301 PetscFunctionBegin; 302 ierr = PetscNewLog(pc,PC_Cholesky,&dir);CHKERRQ(ierr); 303 304 ((PC_Factor*)dir)->fact = 0; 305 dir->inplace = PETSC_FALSE; 306 ierr = MatFactorInfoInitialize(&((PC_Factor*)dir)->info);CHKERRQ(ierr); 307 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 308 ((PC_Factor*)dir)->info.fill = 5.0; 309 ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE; 310 ((PC_Factor*)dir)->info.shiftamount = 0.0; 311 ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON; 312 ((PC_Factor*)dir)->info.pivotinblocks = 1.0; 313 dir->col = 0; 314 dir->row = 0; 315 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 316 ierr = PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 317 dir->reusefill = PETSC_FALSE; 318 dir->reuseordering = PETSC_FALSE; 319 pc->data = (void*)dir; 320 321 pc->ops->destroy = PCDestroy_Cholesky; 322 pc->ops->reset = PCReset_Cholesky; 323 pc->ops->apply = PCApply_Cholesky; 324 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 325 pc->ops->setup = PCSetUp_Cholesky; 326 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 327 pc->ops->view = PCView_Cholesky; 328 pc->ops->applyrichardson = 0; 329 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 330 331 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor", 332 PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr); 333 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 334 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 335 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 336 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 337 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 338 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 339 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor", 340 PCFactorSetShiftType_Factor);CHKERRQ(ierr); 341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor", 342 PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 343 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 344 PCFactorSetFill_Factor);CHKERRQ(ierr); 345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky", 346 PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr); 347 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 348 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 349 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky", 350 PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr); 351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky", 352 PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr); 353 PetscFunctionReturn(0); 354 } 355 EXTERN_C_END 356