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 IS row,col; /* index sets used for reordering */ 12 } PC_Cholesky; 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "PCSetFromOptions_Cholesky" 16 static PetscErrorCode PCSetFromOptions_Cholesky(PetscOptionItems *PetscOptionsObject,PC pc) 17 { 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 ierr = PetscOptionsHead(PetscOptionsObject,"Cholesky options");CHKERRQ(ierr); 22 ierr = PCSetFromOptions_Factor(PetscOptionsObject,pc);CHKERRQ(ierr); 23 ierr = PetscOptionsTail();CHKERRQ(ierr); 24 PetscFunctionReturn(0); 25 } 26 27 #undef __FUNCT__ 28 #define __FUNCT__ "PCView_Cholesky" 29 static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer) 30 { 31 PC_Cholesky *chol = (PC_Cholesky*)pc->data; 32 PetscErrorCode ierr; 33 PetscBool iascii; 34 35 PetscFunctionBegin; 36 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 37 if (iascii) { 38 if (chol->hdr.inplace) { 39 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");CHKERRQ(ierr); 40 } else { 41 ierr = PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");CHKERRQ(ierr); 42 } 43 44 if (chol->hdr.reusefill) {ierr = PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");CHKERRQ(ierr);} 45 if (chol->hdr.reuseordering) {ierr = PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");CHKERRQ(ierr);} 46 } 47 ierr = PCView_Factor(pc,viewer);CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "PCSetUp_Cholesky" 53 static PetscErrorCode PCSetUp_Cholesky(PC pc) 54 { 55 PetscErrorCode ierr; 56 PetscBool flg; 57 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 58 const MatSolverPackage stype; 59 MatFactorError err; 60 61 PetscFunctionBegin; 62 pc->failedreason = PC_NOERROR; 63 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill; 64 65 ierr = MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);CHKERRQ(ierr); 66 if (dir->hdr.inplace) { 67 if (dir->row && dir->col && (dir->row != dir->col)) { 68 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 69 } 70 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 71 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 72 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */ 73 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 74 } 75 if (dir->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr);} 76 ierr = MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 77 ierr = MatFactorGetError(pc->pmat,&err);CHKERRQ(ierr); 78 if (err) { /* Factor() fails */ 79 pc->failedreason = (PCFailedReason)err; 80 PetscFunctionReturn(0); 81 } 82 83 ((PC_Factor*)dir)->fact = pc->pmat; 84 } else { 85 MatInfo info; 86 87 if (!pc->setupcalled) { 88 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 89 /* check if dir->row == dir->col */ 90 ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr); 91 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal"); 92 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */ 93 94 flg = PETSC_FALSE; 95 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);CHKERRQ(ierr); 96 if (flg) { 97 PetscReal tol = 1.e-10; 98 ierr = PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);CHKERRQ(ierr); 99 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 100 } 101 if (dir->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr);} 102 if (!((PC_Factor*)dir)->fact) { 103 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 104 } 105 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 106 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 107 dir->hdr.actualfill = info.fill_ratio_needed; 108 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);CHKERRQ(ierr); 109 } else if (pc->flag != SAME_NONZERO_PATTERN) { 110 if (!dir->hdr.reuseordering) { 111 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 112 ierr = MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);CHKERRQ(ierr); 113 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); /* only use dir->row ordering in CholeskyFactor */ 114 115 flg = PETSC_FALSE; 116 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);CHKERRQ(ierr); 117 if (flg) { 118 PetscReal tol = 1.e-10; 119 ierr = PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);CHKERRQ(ierr); 120 ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr); 121 } 122 if (dir->row) {ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);CHKERRQ(ierr);} 123 } 124 ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 125 ierr = MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);CHKERRQ(ierr); 126 ierr = MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 127 ierr = MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 128 dir->hdr.actualfill = info.fill_ratio_needed; 129 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);CHKERRQ(ierr); 130 } else { 131 ierr = MatFactorGetError(((PC_Factor*)dir)->fact,&err);CHKERRQ(ierr); 132 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) { 133 ierr = MatFactorClearError(((PC_Factor*)dir)->fact);CHKERRQ(ierr); 134 pc->failedreason = PC_NOERROR; 135 } 136 } 137 ierr = MatFactorGetError(((PC_Factor*)dir)->fact,&err);CHKERRQ(ierr); 138 if (err) { /* FactorSymbolic() fails */ 139 pc->failedreason = (PCFailedReason)err; 140 PetscFunctionReturn(0); 141 } 142 143 ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);CHKERRQ(ierr); 144 ierr = MatFactorGetError(((PC_Factor*)dir)->fact,&err);CHKERRQ(ierr); 145 if (err) { /* FactorNumeric() fails */ 146 pc->failedreason = (PCFailedReason)err; 147 } 148 } 149 150 ierr = PCFactorGetMatSolverPackage(pc,&stype);CHKERRQ(ierr); 151 if (!stype) { 152 const MatSolverPackage solverpackage; 153 ierr = MatFactorGetSolverPackage(((PC_Factor*)dir)->fact,&solverpackage);CHKERRQ(ierr); 154 ierr = PCFactorSetMatSolverPackage(pc,solverpackage);CHKERRQ(ierr); 155 } 156 PetscFunctionReturn(0); 157 } 158 159 #undef __FUNCT__ 160 #define __FUNCT__ "PCReset_Cholesky" 161 static PetscErrorCode PCReset_Cholesky(PC pc) 162 { 163 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {ierr = MatDestroy(&((PC_Factor*)dir)->fact);CHKERRQ(ierr);} 168 ierr = ISDestroy(&dir->row);CHKERRQ(ierr); 169 ierr = ISDestroy(&dir->col);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "PCDestroy_Cholesky" 175 static PetscErrorCode PCDestroy_Cholesky(PC pc) 176 { 177 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 ierr = PCReset_Cholesky(pc);CHKERRQ(ierr); 182 ierr = PetscFree(((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 183 ierr = PetscFree(((PC_Factor*)dir)->solvertype);CHKERRQ(ierr); 184 ierr = PetscFree(pc->data);CHKERRQ(ierr); 185 PetscFunctionReturn(0); 186 } 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "PCApply_Cholesky" 190 static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y) 191 { 192 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 if (dir->hdr.inplace) { 197 ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr); 198 } else { 199 ierr = MatSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 200 } 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->hdr.inplace) { 213 ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr); 214 } else { 215 ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 216 } 217 PetscFunctionReturn(0); 218 } 219 220 /* -----------------------------------------------------------------------------------*/ 221 222 /* -----------------------------------------------------------------------------------*/ 223 224 #undef __FUNCT__ 225 #define __FUNCT__ "PCFactorSetReuseOrdering" 226 /*@ 227 PCFactorSetReuseOrdering - When similar matrices are factored, this 228 causes the ordering computed in the first factor to be used for all 229 following factors. 230 231 Logically Collective on PC 232 233 Input Parameters: 234 + pc - the preconditioner context 235 - flag - PETSC_TRUE to reuse else PETSC_FALSE 236 237 Options Database Key: 238 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 239 240 Level: intermediate 241 242 .keywords: PC, levels, reordering, factorization, incomplete, LU 243 244 .seealso: PCFactorSetReuseFill() 245 @*/ 246 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 247 { 248 PetscErrorCode ierr; 249 250 PetscFunctionBegin; 251 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 252 PetscValidLogicalCollectiveBool(pc,flag,2); 253 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 /*MC 258 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 259 260 Options Database Keys: 261 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 262 . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu 263 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 264 . -pc_factor_fill <fill> - Sets fill amount 265 . -pc_factor_in_place - Activates in-place factorization 266 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 267 268 Notes: Not all options work for all matrix formats 269 270 Level: beginner 271 272 Concepts: Cholesky factorization, direct solver 273 274 Notes: Usually this will compute an "exact" solution in one iteration and does 275 not need a Krylov method (i.e. you can use -ksp_type preonly, or 276 KSPSetType(ksp,KSPPREONLY) for the Krylov method 277 278 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 279 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 280 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 281 PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType() 282 283 M*/ 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "PCCreate_Cholesky" 287 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) 288 { 289 PetscErrorCode ierr; 290 PC_Cholesky *dir; 291 292 PetscFunctionBegin; 293 ierr = PetscNewLog(pc,&dir);CHKERRQ(ierr); 294 pc->data = (void*)dir; 295 ierr = PCFactorInitialize(pc);CHKERRQ(ierr); 296 297 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 298 ((PC_Factor*)dir)->info.fill = 5.0; 299 300 dir->col = 0; 301 dir->row = 0; 302 303 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 304 305 pc->ops->destroy = PCDestroy_Cholesky; 306 pc->ops->reset = PCReset_Cholesky; 307 pc->ops->apply = PCApply_Cholesky; 308 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 309 pc->ops->setup = PCSetUp_Cholesky; 310 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 311 pc->ops->view = PCView_Cholesky; 312 pc->ops->applyrichardson = 0; 313 PetscFunctionReturn(0); 314 } 315