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:   PetscFunctionBegin;
 17:   PetscOptionsHeadBegin(PetscOptionsObject, "Cholesky options");
 18:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
 19:   PetscOptionsHeadEnd();
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

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

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

 35:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
 36:   PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));

 38:   PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
 39:   if (dir->hdr.inplace) {
 40:     MatFactorType ftype;

 42:     PetscCall(MatGetFactorType(pc->pmat, &ftype));
 43:     if (ftype == MAT_FACTOR_NONE) {
 44:       if (dir->row && dir->col && (dir->row != dir->col)) PetscCall(ISDestroy(&dir->row));
 45:       PetscCall(ISDestroy(&dir->col));
 46:       /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
 47:       PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
 48:       PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
 49:       if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
 50:         PetscCall(ISDestroy(&dir->col));
 51:       }
 52:       PetscCall(MatCholeskyFactor(pc->pmat, dir->row, &((PC_Factor *)dir)->info));
 53:       PetscCall(MatFactorGetError(pc->pmat, &err));
 54:       if (err) { /* Factor() fails */
 55:         pc->failedreason = (PCFailedReason)err;
 56:         PetscFunctionReturn(PETSC_SUCCESS);
 57:       }
 58:     }
 59:     ((PC_Factor *)dir)->fact = pc->pmat;
 60:   } else {
 61:     MatInfo info;

 63:     if (!pc->setupcalled) {
 64:       PetscBool canuseordering;
 65:       if (!((PC_Factor *)dir)->fact) { PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_CHOLESKY, &((PC_Factor *)dir)->fact)); }
 66:       PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
 67:       if (canuseordering) {
 68:         PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
 69:         PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
 70:         /* check if dir->row == dir->col */
 71:         if (dir->row) {
 72:           PetscCall(ISEqual(dir->row, dir->col, &flg));
 73:           PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "row and column permutations must be equal");
 74:         }
 75:         PetscCall(ISDestroy(&dir->col)); /* only pass one ordering into CholeskyFactor */

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

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

125:     PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
126:     PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
127:     if (err) { /* FactorNumeric() fails */
128:       pc->failedreason = (PCFailedReason)err;
129:     }
130:   }

132:   PetscCall(PCFactorGetMatSolverType(pc, &stype));
133:   if (!stype) {
134:     MatSolverType solverpackage;
135:     PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
136:     PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
137:   }
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: static PetscErrorCode PCReset_Cholesky(PC pc)
142: {
143:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

145:   PetscFunctionBegin;
146:   if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
147:   PetscCall(ISDestroy(&dir->row));
148:   PetscCall(ISDestroy(&dir->col));
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

152: static PetscErrorCode PCDestroy_Cholesky(PC pc)
153: {
154:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

156:   PetscFunctionBegin;
157:   PetscCall(PCReset_Cholesky(pc));
158:   PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
159:   PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
160:   PetscCall(PCFactorClearComposedFunctions(pc));
161:   PetscCall(PetscFree(pc->data));
162:   PetscFunctionReturn(PETSC_SUCCESS);
163: }

165: static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y)
166: {
167:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

169:   PetscFunctionBegin;
170:   if (dir->hdr.inplace) {
171:     PetscCall(MatSolve(pc->pmat, x, y));
172:   } else {
173:     PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y));
174:   }
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y)
179: {
180:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

182:   PetscFunctionBegin;
183:   if (dir->hdr.inplace) {
184:     PetscCall(MatMatSolve(pc->pmat, X, Y));
185:   } else {
186:     PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y));
187:   }
188:   PetscFunctionReturn(PETSC_SUCCESS);
189: }

191: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y)
192: {
193:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

195:   PetscFunctionBegin;
196:   if (dir->hdr.inplace) {
197:     PetscCall(MatForwardSolve(pc->pmat, x, y));
198:   } else {
199:     PetscCall(MatForwardSolve(((PC_Factor *)dir)->fact, x, y));
200:   }
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

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

208:   PetscFunctionBegin;
209:   if (dir->hdr.inplace) {
210:     PetscCall(MatBackwardSolve(pc->pmat, x, y));
211:   } else {
212:     PetscCall(MatBackwardSolve(((PC_Factor *)dir)->fact, x, y));
213:   }
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y)
218: {
219:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

221:   PetscFunctionBegin;
222:   if (dir->hdr.inplace) {
223:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
224:   } else {
225:     PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y));
226:   }
227:   PetscFunctionReturn(PETSC_SUCCESS);
228: }

230: /*@
231:    PCFactorSetReuseOrdering - When similar matrices are factored, this
232:    causes the ordering computed in the first factor to be used for all
233:    following factors.

235:    Logically Collective

237:    Input Parameters:
238: +  pc - the preconditioner context
239: -  flag - `PETSC_TRUE` to reuse else `PETSC_FALSE`

241:    Options Database Key:
242: .  -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`

244:    Level: intermediate

246: .seealso: `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()`
247: @*/
248: PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag)
249: {
250:   PetscFunctionBegin;
253:   PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

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

260:    Options Database Keys:
261: +  -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
262: .  -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
263: .  -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
264: .  -pc_factor_fill <fill> - Sets fill amount
265: .  -pc_factor_in_place - Activates in-place factorization
266: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

268:    Level: beginner

270:    Notes:
271:    Not all options work for all matrix formats

273:    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

277: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
278:           `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
279:           `PCFactorSetFill()`, `PCFactorSetShiftNonzero()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
280:           `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`
281: M*/

283: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
284: {
285:   PC_Cholesky *dir;

287:   PetscFunctionBegin;
288:   PetscCall(PetscNew(&dir));
289:   pc->data = (void *)dir;
290:   PetscCall(PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY));

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

294:   pc->ops->destroy             = PCDestroy_Cholesky;
295:   pc->ops->reset               = PCReset_Cholesky;
296:   pc->ops->apply               = PCApply_Cholesky;
297:   pc->ops->matapply            = PCMatApply_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     = NULL;
305:   PetscFunctionReturn(PETSC_SUCCESS);
306: }