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