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