Actual source code: cholesky.c


  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>

  9: typedef struct {
 10:   PC_Factor hdr;
 11:   IS        row, col; /* index sets used for reordering */
 12: } PC_Cholesky;

 14: static PetscErrorCode PCSetFromOptions_Cholesky(PC pc, PetscOptionItems *PetscOptionsObject)
 15: {
 16:   PetscOptionsHeadBegin(PetscOptionsObject, "Cholesky options");
 17:   PCSetFromOptions_Factor(pc, PetscOptionsObject);
 18:   PetscOptionsHeadEnd();
 19:   return 0;
 20: }

 22: static PetscErrorCode PCSetUp_Cholesky(PC pc)
 23: {
 24:   PetscBool      flg;
 25:   PC_Cholesky   *dir = (PC_Cholesky *)pc->data;
 26:   MatSolverType  stype;
 27:   MatFactorError err;
 28:   const char    *prefix;

 30:   pc->failedreason = PC_NOERROR;
 31:   if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;

 33:   PCGetOptionsPrefix(pc, &prefix);
 34:   MatSetOptionsPrefixFactor(pc->pmat, prefix);

 36:   MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
 37:   if (dir->hdr.inplace) {
 38:     if (dir->row && dir->col && (dir->row != dir->col)) ISDestroy(&dir->row);
 39:     ISDestroy(&dir->col);
 40:     /* should only get reordering if the factor matrix uses it but cannot determine because MatGetFactor() not called */
 41:     PCFactorSetDefaultOrdering_Factor(pc);
 42:     MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
 43:     if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
 44:       ISDestroy(&dir->col);
 45:     }
 46:     MatCholeskyFactor(pc->pmat, dir->row, &((PC_Factor *)dir)->info);
 47:     MatFactorGetError(pc->pmat, &err);
 48:     if (err) { /* Factor() fails */
 49:       pc->failedreason = (PCFailedReason)err;
 50:       return 0;
 51:     }

 53:     ((PC_Factor *)dir)->fact = pc->pmat;
 54:   } else {
 55:     MatInfo info;

 57:     if (!pc->setupcalled) {
 58:       PetscBool canuseordering;
 59:       if (!((PC_Factor *)dir)->fact) { MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_CHOLESKY, &((PC_Factor *)dir)->fact); }
 60:       MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering);
 61:       if (canuseordering) {
 62:         PCFactorSetDefaultOrdering_Factor(pc);
 63:         MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
 64:         /* check if dir->row == dir->col */
 65:         if (dir->row) {
 66:           ISEqual(dir->row, dir->col, &flg);
 68:         }
 69:         ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */

 71:         flg = PETSC_FALSE;
 72:         PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg, NULL);
 73:         if (flg) {
 74:           PetscReal tol = 1.e-10;
 75:           PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL);
 76:           MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row);
 77:         }
 78:       }
 79:       MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info);
 80:       MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
 81:       dir->hdr.actualfill = info.fill_ratio_needed;
 82:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
 83:       if (!dir->hdr.reuseordering) {
 84:         PetscBool canuseordering;
 85:         MatDestroy(&((PC_Factor *)dir)->fact);
 86:         MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_CHOLESKY, &((PC_Factor *)dir)->fact);
 87:         MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering);
 88:         if (canuseordering) {
 89:           ISDestroy(&dir->row);
 90:           PCFactorSetDefaultOrdering_Factor(pc);
 91:           MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
 92:           ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */

 94:           flg = PETSC_FALSE;
 95:           PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg, NULL);
 96:           if (flg) {
 97:             PetscReal tol = 1.e-10;
 98:             PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL);
 99:             MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row);
100:           }
101:         }
102:       }
103:       MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info);
104:       MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
105:       dir->hdr.actualfill = info.fill_ratio_needed;
106:     } else {
107:       MatFactorGetError(((PC_Factor *)dir)->fact, &err);
108:       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
109:         MatFactorClearError(((PC_Factor *)dir)->fact);
110:         pc->failedreason = PC_NOERROR;
111:       }
112:     }
113:     MatFactorGetError(((PC_Factor *)dir)->fact, &err);
114:     if (err) { /* FactorSymbolic() fails */
115:       pc->failedreason = (PCFailedReason)err;
116:       return 0;
117:     }

119:     MatCholeskyFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info);
120:     MatFactorGetError(((PC_Factor *)dir)->fact, &err);
121:     if (err) { /* FactorNumeric() fails */
122:       pc->failedreason = (PCFailedReason)err;
123:     }
124:   }

126:   PCFactorGetMatSolverType(pc, &stype);
127:   if (!stype) {
128:     MatSolverType solverpackage;
129:     MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage);
130:     PCFactorSetMatSolverType(pc, solverpackage);
131:   }
132:   return 0;
133: }

135: static PetscErrorCode PCReset_Cholesky(PC pc)
136: {
137:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

139:   if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) MatDestroy(&((PC_Factor *)dir)->fact);
140:   ISDestroy(&dir->row);
141:   ISDestroy(&dir->col);
142:   return 0;
143: }

145: static PetscErrorCode PCDestroy_Cholesky(PC pc)
146: {
147:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

149:   PCReset_Cholesky(pc);
150:   PetscFree(((PC_Factor *)dir)->ordering);
151:   PetscFree(((PC_Factor *)dir)->solvertype);
152:   PCFactorClearComposedFunctions(pc);
153:   PetscFree(pc->data);
154:   return 0;
155: }

157: static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y)
158: {
159:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

161:   if (dir->hdr.inplace) {
162:     MatSolve(pc->pmat, x, y);
163:   } else {
164:     MatSolve(((PC_Factor *)dir)->fact, x, y);
165:   }
166:   return 0;
167: }

169: static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y)
170: {
171:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

173:   if (dir->hdr.inplace) {
174:     MatMatSolve(pc->pmat, X, Y);
175:   } else {
176:     MatMatSolve(((PC_Factor *)dir)->fact, X, Y);
177:   }
178:   return 0;
179: }

181: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y)
182: {
183:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

185:   if (dir->hdr.inplace) {
186:     MatForwardSolve(pc->pmat, x, y);
187:   } else {
188:     MatForwardSolve(((PC_Factor *)dir)->fact, x, y);
189:   }
190:   return 0;
191: }

193: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc, Vec x, Vec y)
194: {
195:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

197:   if (dir->hdr.inplace) {
198:     MatBackwardSolve(pc->pmat, x, y);
199:   } else {
200:     MatBackwardSolve(((PC_Factor *)dir)->fact, x, y);
201:   }
202:   return 0;
203: }

205: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y)
206: {
207:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

209:   if (dir->hdr.inplace) {
210:     MatSolveTranspose(pc->pmat, x, y);
211:   } else {
212:     MatSolveTranspose(((PC_Factor *)dir)->fact, x, y);
213:   }
214:   return 0;
215: }

217: /*@
218:    PCFactorSetReuseOrdering - When similar matrices are factored, this
219:    causes the ordering computed in the first factor to be used for all
220:    following factors.

222:    Logically Collective

224:    Input Parameters:
225: +  pc - the preconditioner context
226: -  flag - `PETSC_TRUE` to reuse else `PETSC_FALSE`

228:    Options Database Key:
229: .  -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`

231:    Level: intermediate

233: .seealso: `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()`
234: @*/
235: PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag)
236: {
239:   PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag));
240:   return 0;
241: }

243: /*MC
244:    PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner

246:    Options Database Keys:
247: +  -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
248: .  -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
249: .  -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
250: .  -pc_factor_fill <fill> - Sets fill amount
251: .  -pc_factor_in_place - Activates in-place factorization
252: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

254:    Level: beginner

256:    Notes:
257:    Not all options work for all matrix formats

259:    Usually this will compute an "exact" solution in one iteration and does
260:    not need a Krylov method (i.e. you can use -ksp_type preonly, or
261:    `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method

263: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
264:           `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
265:           `PCFactorSetFill()`, `PCFactorSetShiftNonzero()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
266:           `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`
267: M*/

269: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
270: {
271:   PC_Cholesky *dir;

273:   PetscNew(&dir);
274:   pc->data = (void *)dir;
275:   PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY);

277:   ((PC_Factor *)dir)->info.fill = 5.0;

279:   pc->ops->destroy             = PCDestroy_Cholesky;
280:   pc->ops->reset               = PCReset_Cholesky;
281:   pc->ops->apply               = PCApply_Cholesky;
282:   pc->ops->matapply            = PCMatApply_Cholesky;
283:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Cholesky;
284:   pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
285:   pc->ops->applytranspose      = PCApplyTranspose_Cholesky;
286:   pc->ops->setup               = PCSetUp_Cholesky;
287:   pc->ops->setfromoptions      = PCSetFromOptions_Cholesky;
288:   pc->ops->view                = PCView_Factor;
289:   pc->ops->applyrichardson     = NULL;
290:   return 0;
291: }