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