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