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