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 107 ((PC_Factor*)dir)->fact = pc->pmat; 108 } else { 109 MatInfo info; 110 if (!pc->setupcalled) { 111 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 112 /* check if dir->row == dir->col */ 113 ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr); 114 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal"); 115 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */ 116 117 flg = PETSC_FALSE; 118 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 119 if (flg) { 120 PetscReal tol = 1.e-10; 121 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 122 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 123 } 124 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 125 if (!((PC_Factor*)dir)->fact) { 126 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 127 } 128 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 129 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 130 dir->actualfill = info.fill_ratio_needed; 131 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 132 } else if (pc->flag != SAME_NONZERO_PATTERN) { 133 if (!dir->reuseordering) { 134 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 135 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 136 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only use dir->row ordering in CholeskyFactor */ 137 138 flg = PETSC_FALSE; 139 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 140 if (flg) { 141 PetscReal tol = 1.e-10; 142 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 143 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 144 } 145 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 146 } 147 ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 148 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 149 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 150 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 151 dir->actualfill = info.fill_ratio_needed; 152 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 153 } 154 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 155 } 156 PetscFunctionReturn(0); 157 } 158 159 #undef __FUNCT__ 160 #define __FUNCT__ "PCReset_Cholesky" 161 static PetscErrorCode PCReset_Cholesky(PC pc) 162 { 163 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 if (!dir->inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 168 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 169 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "PCDestroy_Cholesky" 175 static PetscErrorCode PCDestroy_Cholesky(PC pc) 176 { 177 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 ierr = PCReset_Cholesky(pc);CHKERRQ(ierr); 182 ierr = PetscFree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 183 ierr = PetscFree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 184 ierr = PetscFree(pc->data);CHKERRQ(ierr); 185 PetscFunctionReturn(0); 186 } 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "PCApply_Cholesky" 190 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 191 { 192 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 if (dir->inplace) { 197 ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr); 198 } else { 199 ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 200 } 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) { 213 ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr); 214 } else { 215 ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 216 } 217 PetscFunctionReturn(0); 218 } 219 220 /* -----------------------------------------------------------------------------------*/ 221 222 EXTERN_C_BEGIN 223 #undef __FUNCT__ 224 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky" 225 PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc) 226 { 227 PC_Cholesky *dir; 228 229 PetscFunctionBegin; 230 dir = (PC_Cholesky*)pc->data; 231 dir->inplace = PETSC_TRUE; 232 PetscFunctionReturn(0); 233 } 234 EXTERN_C_END 235 236 /* -----------------------------------------------------------------------------------*/ 237 238 #undef __FUNCT__ 239 #define __FUNCT__ "PCFactorSetReuseOrdering" 240 /*@ 241 PCFactorSetReuseOrdering - When similar matrices are factored, this 242 causes the ordering computed in the first factor to be used for all 243 following factors. 244 245 Logically Collective on PC 246 247 Input Parameters: 248 + pc - the preconditioner context 249 - flag - PETSC_TRUE to reuse else PETSC_FALSE 250 251 Options Database Key: 252 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 253 254 Level: intermediate 255 256 .keywords: PC, levels, reordering, factorization, incomplete, LU 257 258 .seealso: PCFactorSetReuseFill() 259 @*/ 260 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 261 { 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 266 PetscValidLogicalCollectiveBool(pc,flag,2); 267 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 271 272 /*MC 273 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 274 275 Options Database Keys: 276 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 277 . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu 278 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 279 . -pc_factor_fill <fill> - Sets fill amount 280 . -pc_factor_in_place - Activates in-place factorization 281 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 282 283 Notes: Not all options work for all matrix formats 284 285 Level: beginner 286 287 Concepts: Cholesky factorization, direct solver 288 289 Notes: Usually this will compute an "exact" solution in one iteration and does 290 not need a Krylov method (i.e. you can use -ksp_type preonly, or 291 KSPSetType(ksp,KSPPREONLY) for the Krylov method 292 293 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 294 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 295 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 296 PCFactorSetUseInPlace(), PCFactorSetMatOrderingType() 297 298 M*/ 299 300 EXTERN_C_BEGIN 301 #undef __FUNCT__ 302 #define __FUNCT__ "PCCreate_Cholesky" 303 PetscErrorCode PCCreate_Cholesky(PC pc) 304 { 305 PetscErrorCode ierr; 306 PC_Cholesky *dir; 307 308 PetscFunctionBegin; 309 ierr = PetscNewLog(pc,PC_Cholesky,&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 = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor", 344 PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr); 345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 346 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 347 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 348 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 349 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 350 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor", 352 PCFactorSetShiftType_Factor);CHKERRQ(ierr); 353 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor", 354 PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 356 PCFactorSetFill_Factor);CHKERRQ(ierr); 357 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky", 358 PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr); 359 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 360 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 361 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky", 362 PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr); 363 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky", 364 PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367 EXTERN_C_END 368