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