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