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 PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y) 171 { 172 PC_Cholesky *dir = (PC_Cholesky*)pc->data; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 if (dir->hdr.inplace) { 177 ierr = MatForwardSolve(pc->pmat,x,y);CHKERRQ(ierr); 178 } else { 179 ierr = MatForwardSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 180 } 181 PetscFunctionReturn(0); 182 } 183 184 static PetscErrorCode PCApplySymmetricRight_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 = MatBackwardSolve(pc->pmat,x,y);CHKERRQ(ierr); 192 } else { 193 ierr = MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 static PetscErrorCode PCApplyTranspose_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 = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr); 206 } else { 207 ierr = MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);CHKERRQ(ierr); 208 } 209 PetscFunctionReturn(0); 210 } 211 212 /* -----------------------------------------------------------------------------------*/ 213 214 /* -----------------------------------------------------------------------------------*/ 215 216 /*@ 217 PCFactorSetReuseOrdering - When similar matrices are factored, this 218 causes the ordering computed in the first factor to be used for all 219 following factors. 220 221 Logically Collective on PC 222 223 Input Parameters: 224 + pc - the preconditioner context 225 - flag - PETSC_TRUE to reuse else PETSC_FALSE 226 227 Options Database Key: 228 . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 229 230 Level: intermediate 231 232 .seealso: PCFactorSetReuseFill() 233 @*/ 234 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 235 { 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 240 PetscValidLogicalCollectiveBool(pc,flag,2); 241 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244 245 /*MC 246 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 247 248 Options Database Keys: 249 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 250 . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu 251 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 252 . -pc_factor_fill <fill> - Sets fill amount 253 . -pc_factor_in_place - Activates in-place factorization 254 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 255 256 Notes: 257 Not all options work for all matrix formats 258 259 Level: beginner 260 261 Notes: 262 Usually this will compute an "exact" solution in one iteration and does 263 not need a Krylov method (i.e. you can use -ksp_type preonly, or 264 KSPSetType(ksp,KSPPREONLY) for the Krylov method 265 266 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 267 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 268 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 269 PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType() 270 271 M*/ 272 273 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) 274 { 275 PetscErrorCode ierr; 276 PC_Cholesky *dir; 277 278 PetscFunctionBegin; 279 ierr = PetscNewLog(pc,&dir);CHKERRQ(ierr); 280 pc->data = (void*)dir; 281 ierr = PCFactorInitialize(pc);CHKERRQ(ierr); 282 283 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 284 ((PC_Factor*)dir)->info.fill = 5.0; 285 286 dir->col = NULL; 287 dir->row = NULL; 288 289 /* MATORDERINGNATURAL_OR_ND allows selecting type based on matrix type sbaij or aij */ 290 ierr = PetscStrallocpy(MATORDERINGNATURAL_OR_ND,(char**)&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 291 292 pc->ops->destroy = PCDestroy_Cholesky; 293 pc->ops->reset = PCReset_Cholesky; 294 pc->ops->apply = PCApply_Cholesky; 295 pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky; 296 pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky; 297 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 298 pc->ops->setup = PCSetUp_Cholesky; 299 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 300 pc->ops->view = PCView_Factor; 301 pc->ops->applyrichardson = NULL; 302 PetscFunctionReturn(0); 303 } 304