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