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