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