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,PETSC_VIEWER_ASCII,&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_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 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 132 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 133 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 134 dir->actualfill = info.fill_ratio_needed; 135 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 136 } else if (pc->flag != SAME_NONZERO_PATTERN) { 137 if (!dir->reuseordering) { 138 if (dir->row && dir->col && (dir->row != dir->col)) { 139 ierr = ISDestroy(dir->row);CHKERRQ(ierr); 140 dir->row = 0; 141 } 142 if (dir->col) { 143 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 144 dir->col =0; 145 } 146 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 147 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 148 ierr = ISDestroy(dir->col);CHKERRQ(ierr); 149 dir->col=0; 150 } 151 flg = PETSC_FALSE; 152 ierr = PetscOptionsGetTruth(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);CHKERRQ(ierr); 153 if (flg) { 154 PetscReal tol = 1.e-10; 155 ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr); 156 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 157 } 158 if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);} 159 } 160 ierr = MatDestroy(((PC_Factor*)dir)->fact);CHKERRQ(ierr); 161 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 162 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 163 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 164 dir->actualfill = info.fill_ratio_needed; 165 ierr = PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);CHKERRQ(ierr); 166 } 167 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 168 } 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCDestroy_Cholesky" 174 static PetscErrorCode PCDestroy_Cholesky(PC pc) 175 { 176 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 177 PetscErrorCode ierr; 178 179 PetscFunctionBegin; 180 if (!dir->inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 181 if (dir->row) {ierr = ISDestroy(dir->row);CHKERRQ(ierr);} 182 if (dir->col) {ierr = ISDestroy(dir->col);CHKERRQ(ierr);} 183 ierr = PetscFree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 184 ierr = PetscFree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 185 ierr = PetscFree(dir);CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 #undef __FUNCT__ 190 #define __FUNCT__ "PCApply_Cholesky" 191 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 192 { 193 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 if (dir->inplace) {ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr);} 198 else {ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "PCApplyTranspose_Cholesky" 204 static PetscErrorCode PCApplyTranspose_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) {ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr);} 211 else {ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr);} 212 PetscFunctionReturn(0); 213 } 214 215 /* -----------------------------------------------------------------------------------*/ 216 217 EXTERN_C_BEGIN 218 #undef __FUNCT__ 219 #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky" 220 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseInPlace_Cholesky(PC pc) 221 { 222 PC_Cholesky *dir; 223 224 PetscFunctionBegin; 225 dir = (PC_Cholesky*)pc->data; 226 dir->inplace = PETSC_TRUE; 227 PetscFunctionReturn(0); 228 } 229 EXTERN_C_END 230 231 /* -----------------------------------------------------------------------------------*/ 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "PCFactorSetReuseOrdering" 235 /*@ 236 PCFactorSetReuseOrdering - When similar matrices are factored, this 237 causes the ordering computed in the first factor to be used for all 238 following factors. 239 240 Collective on PC 241 242 Input Parameters: 243 + pc - the preconditioner context 244 - flag - PETSC_TRUE to reuse else PETSC_FALSE 245 246 Options Database Key: 247 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 248 249 Level: intermediate 250 251 .keywords: PC, levels, reordering, factorization, incomplete, LU 252 253 .seealso: PCFactorSetReuseFill() 254 @*/ 255 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering(PC pc,PetscTruth flag) 256 { 257 PetscErrorCode ierr,(*f)(PC,PetscTruth); 258 259 PetscFunctionBegin; 260 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 261 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",(void (**)(void))&f);CHKERRQ(ierr); 262 if (f) { 263 ierr = (*f)(pc,flag);CHKERRQ(ierr); 264 } 265 PetscFunctionReturn(0); 266 } 267 268 269 /*MC 270 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 271 272 Options Database Keys: 273 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 274 . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles 275 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 276 . -pc_factor_fill <fill> - Sets fill amount 277 . -pc_factor_in_place - Activates in-place factorization 278 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 279 280 Notes: Not all options work for all matrix formats 281 282 Level: beginner 283 284 Concepts: Cholesky factorization, direct solver 285 286 Notes: Usually this will compute an "exact" solution in one iteration and does 287 not need a Krylov method (i.e. you can use -ksp_type preonly, or 288 KSPSetType(ksp,KSPPREONLY) for the Krylov method 289 290 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 291 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 292 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 293 PCFactorSetUseInPlace(), PCFactorSetMatOrderingType() 294 295 M*/ 296 297 EXTERN_C_BEGIN 298 #undef __FUNCT__ 299 #define __FUNCT__ "PCCreate_Cholesky" 300 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Cholesky(PC pc) 301 { 302 PetscErrorCode ierr; 303 PC_Cholesky *dir; 304 305 PetscFunctionBegin; 306 ierr = PetscNewLog(pc,PC_Cholesky,&dir);CHKERRQ(ierr); 307 308 ((PC_Factor*)dir)->fact = 0; 309 dir->inplace = PETSC_FALSE; 310 ierr = MatFactorInfoInitialize(&((PC_Factor*)dir)->info);CHKERRQ(ierr); 311 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 312 ((PC_Factor*)dir)->info.fill = 5.0; 313 ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE; 314 ((PC_Factor*)dir)->info.shiftamount = 0.0; 315 ((PC_Factor*)dir)->info.zeropivot = 1.e-12; 316 ((PC_Factor*)dir)->info.pivotinblocks = 1.0; 317 dir->col = 0; 318 dir->row = 0; 319 ierr = PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 320 ierr = PetscStrallocpy(MAT_SOLVER_PETSC,&((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 321 dir->reusefill = PETSC_FALSE; 322 dir->reuseordering = PETSC_FALSE; 323 pc->data = (void*)dir; 324 325 pc->ops->destroy = PCDestroy_Cholesky; 326 pc->ops->apply = PCApply_Cholesky; 327 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 328 pc->ops->setup = PCSetUp_Cholesky; 329 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 330 pc->ops->view = PCView_Cholesky; 331 pc->ops->applyrichardson = 0; 332 pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor; 333 334 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor", 335 PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr); 336 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor", 337 PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr); 338 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor", 339 PCFactorSetZeroPivot_Factor);CHKERRQ(ierr); 340 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor", 341 PCFactorSetShiftType_Factor);CHKERRQ(ierr); 342 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor", 343 PCFactorSetShiftAmount_Factor);CHKERRQ(ierr); 344 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor", 345 PCFactorSetFill_Factor);CHKERRQ(ierr); 346 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky", 347 PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr); 348 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor", 349 PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr); 350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky", 351 PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr); 352 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky", 353 PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr); 354 PetscFunctionReturn(0); 355 } 356 EXTERN_C_END 357