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