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 .keywords: PC, levels, reordering, factorization, incomplete, LU 233 234 .seealso: PCFactorSetReuseFill() 235 @*/ 236 PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag) 237 { 238 PetscErrorCode ierr; 239 240 PetscFunctionBegin; 241 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 242 PetscValidLogicalCollectiveBool(pc,flag,2); 243 ierr = PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246 247 /*MC 248 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner 249 250 Options Database Keys: 251 + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering() 252 . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu 253 . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill() 254 . -pc_factor_fill <fill> - Sets fill amount 255 . -pc_factor_in_place - Activates in-place factorization 256 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine 257 258 Notes: 259 Not all options work for all matrix formats 260 261 Level: beginner 262 263 Concepts: Cholesky factorization, direct solver 264 265 Notes: 266 Usually this will compute an "exact" solution in one iteration and does 267 not need a Krylov method (i.e. you can use -ksp_type preonly, or 268 KSPSetType(ksp,KSPPREONLY) for the Krylov method 269 270 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 271 PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 272 PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount() 273 PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType() 274 275 M*/ 276 277 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc) 278 { 279 PetscErrorCode ierr; 280 PC_Cholesky *dir; 281 282 PetscFunctionBegin; 283 ierr = PetscNewLog(pc,&dir);CHKERRQ(ierr); 284 pc->data = (void*)dir; 285 ierr = PCFactorInitialize(pc);CHKERRQ(ierr); 286 287 ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY; 288 ((PC_Factor*)dir)->info.fill = 5.0; 289 290 dir->col = 0; 291 dir->row = 0; 292 293 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);CHKERRQ(ierr); 294 295 pc->ops->destroy = PCDestroy_Cholesky; 296 pc->ops->reset = PCReset_Cholesky; 297 pc->ops->apply = PCApply_Cholesky; 298 pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky; 299 pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky; 300 pc->ops->applytranspose = PCApplyTranspose_Cholesky; 301 pc->ops->setup = PCSetUp_Cholesky; 302 pc->ops->setfromoptions = PCSetFromOptions_Cholesky; 303 pc->ops->view = PCView_Factor; 304 pc->ops->applyrichardson = 0; 305 PetscFunctionReturn(0); 306 } 307