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