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