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