Actual source code: icc.c
petsc-3.11.4 2019-09-28
2: #include <../src/ksp/pc/impls/factor/icc/icc.h>
4: static PetscErrorCode PCSetUp_ICC(PC pc)
5: {
6: PC_ICC *icc = (PC_ICC*)pc->data;
7: IS perm,cperm;
8: PetscErrorCode ierr;
9: MatInfo info;
10: MatSolverType stype;
11: MatFactorError err;
14: pc->failedreason = PC_NOERROR;
15: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
17: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
18: if (!pc->setupcalled) {
19: if (!((PC_Factor*)icc)->fact) {
20: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
21: }
22: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
23: } else if (pc->flag != SAME_NONZERO_PATTERN) {
24: MatDestroy(&((PC_Factor*)icc)->fact);
25: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
26: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
27: }
28: MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
29: icc->hdr.actualfill = info.fill_ratio_needed;
31: ISDestroy(&cperm);
32: ISDestroy(&perm);
34: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
35: if (err) { /* FactorSymbolic() fails */
36: pc->failedreason = (PCFailedReason)err;
37: return(0);
38: }
39:
40: MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
41: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
42: if (err) { /* FactorNumeric() fails */
43: pc->failedreason = (PCFailedReason)err;
44: }
46: PCFactorGetMatSolverType(pc,&stype);
47: if (!stype) {
48: MatSolverType solverpackage;
49: MatFactorGetSolverType(((PC_Factor*)icc)->fact,&solverpackage);
50: PCFactorSetMatSolverType(pc,solverpackage);
51: }
52: return(0);
53: }
55: static PetscErrorCode PCReset_ICC(PC pc)
56: {
57: PC_ICC *icc = (PC_ICC*)pc->data;
61: MatDestroy(&((PC_Factor*)icc)->fact);
62: return(0);
63: }
65: static PetscErrorCode PCDestroy_ICC(PC pc)
66: {
67: PC_ICC *icc = (PC_ICC*)pc->data;
71: PCReset_ICC(pc);
72: PetscFree(((PC_Factor*)icc)->ordering);
73: PetscFree(((PC_Factor*)icc)->solvertype);
74: PetscFree(pc->data);
75: return(0);
76: }
78: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
79: {
80: PC_ICC *icc = (PC_ICC*)pc->data;
84: MatSolve(((PC_Factor*)icc)->fact,x,y);
85: return(0);
86: }
88: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
89: {
91: PC_ICC *icc = (PC_ICC*)pc->data;
94: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
95: return(0);
96: }
98: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
99: {
101: PC_ICC *icc = (PC_ICC*)pc->data;
104: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
105: return(0);
106: }
108: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
109: {
110: PC_ICC *icc = (PC_ICC*)pc->data;
111: PetscBool flg;
113: /* PetscReal dt[3];*/
116: PetscOptionsHead(PetscOptionsObject,"ICC Options");
117: PCSetFromOptions_Factor(PetscOptionsObject,pc);
119: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
120: /*dt[0] = ((PC_Factor*)icc)->info.dt;
121: dt[1] = ((PC_Factor*)icc)->info.dtcol;
122: dt[2] = ((PC_Factor*)icc)->info.dtcount;
123: PetscInt dtmax = 3;
124: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
125: if (flg) {
126: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
127: }
128: */
129: PetscOptionsTail();
130: return(0);
131: }
133: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
134: {
138: PCView_Factor(pc,viewer);
139: return(0);
140: }
142: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
144: /*MC
145: PCICC - Incomplete Cholesky factorization preconditioners.
147: Options Database Keys:
148: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
149: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
150: its factorization (overwrites original matrix)
151: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
152: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
154: Level: beginner
156: Concepts: incomplete Cholesky factorization
158: Notes:
159: Only implemented for some matrix formats. Not implemented in parallel.
161: For BAIJ matrices this implements a point block ICC.
163: The Manteuffel shift is only implemented for matrices with block size 1
165: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
166: to turn off the shift.
168: References:
169: . 1. - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
170: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
171: Science and Engineering, Kluwer.
174: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
175: PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
176: PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
177: PCFactorSetLevels()
179: M*/
181: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
182: {
184: PC_ICC *icc;
187: PetscNewLog(pc,&icc);
188: pc->data = (void*)icc;
189: PCFactorInitialize(pc);
190: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);
192: ((PC_Factor*)icc)->factortype = MAT_FACTOR_ICC;
193: ((PC_Factor*)icc)->info.fill = 1.0;
194: ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT;
195: ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
197: pc->ops->apply = PCApply_ICC;
198: pc->ops->applytranspose = PCApply_ICC;
199: pc->ops->setup = PCSetUp_ICC;
200: pc->ops->reset = PCReset_ICC;
201: pc->ops->destroy = PCDestroy_ICC;
202: pc->ops->setfromoptions = PCSetFromOptions_ICC;
203: pc->ops->view = PCView_ICC;
204: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
205: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
206: return(0);
207: }