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