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