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