Actual source code: icc.c
petsc-3.14.6 2021-03-30
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 = NULL,cperm = NULL;
8: PetscErrorCode ierr;
9: MatInfo info;
10: MatSolverType stype;
11: MatFactorError err;
14: pc->failedreason = PC_NOERROR;
16: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
17: if (!pc->setupcalled) {
18: if (!((PC_Factor*)icc)->fact) {
19: PetscBool useordering;
20: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
21: MatFactorGetUseOrdering(((PC_Factor*)icc)->fact,&useordering);
22: if (useordering) {
23: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
24: }
25: }
26: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
27: } else if (pc->flag != SAME_NONZERO_PATTERN) {
28: PetscBool useordering;
29: MatDestroy(&((PC_Factor*)icc)->fact);
30: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
31: MatFactorGetUseOrdering(((PC_Factor*)icc)->fact,&useordering);
32: if (useordering) {
33: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
34: }
35: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
36: }
37: MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
38: icc->hdr.actualfill = info.fill_ratio_needed;
40: ISDestroy(&cperm);
41: ISDestroy(&perm);
43: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
44: if (err) { /* FactorSymbolic() fails */
45: pc->failedreason = (PCFailedReason)err;
46: return(0);
47: }
49: MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
50: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
51: if (err) { /* FactorNumeric() fails */
52: pc->failedreason = (PCFailedReason)err;
53: }
55: PCFactorGetMatSolverType(pc,&stype);
56: if (!stype) {
57: MatSolverType solverpackage;
58: MatFactorGetSolverType(((PC_Factor*)icc)->fact,&solverpackage);
59: PCFactorSetMatSolverType(pc,solverpackage);
60: }
61: return(0);
62: }
64: static PetscErrorCode PCReset_ICC(PC pc)
65: {
66: PC_ICC *icc = (PC_ICC*)pc->data;
70: MatDestroy(&((PC_Factor*)icc)->fact);
71: return(0);
72: }
74: static PetscErrorCode PCDestroy_ICC(PC pc)
75: {
76: PC_ICC *icc = (PC_ICC*)pc->data;
80: PCReset_ICC(pc);
81: PetscFree(((PC_Factor*)icc)->ordering);
82: PetscFree(((PC_Factor*)icc)->solvertype);
83: PetscFree(pc->data);
84: return(0);
85: }
87: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
88: {
89: PC_ICC *icc = (PC_ICC*)pc->data;
93: MatSolve(((PC_Factor*)icc)->fact,x,y);
94: return(0);
95: }
97: static PetscErrorCode PCMatApply_ICC(PC pc,Mat X,Mat Y)
98: {
99: PC_ICC *icc = (PC_ICC*)pc->data;
103: MatMatSolve(((PC_Factor*)icc)->fact,X,Y);
104: return(0);
105: }
107: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
108: {
110: PC_ICC *icc = (PC_ICC*)pc->data;
113: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
114: return(0);
115: }
117: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
118: {
120: PC_ICC *icc = (PC_ICC*)pc->data;
123: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
124: return(0);
125: }
127: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
128: {
129: PC_ICC *icc = (PC_ICC*)pc->data;
130: PetscBool flg;
132: /* PetscReal dt[3];*/
135: PetscOptionsHead(PetscOptionsObject,"ICC Options");
136: PCSetFromOptions_Factor(PetscOptionsObject,pc);
138: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
139: /*dt[0] = ((PC_Factor*)icc)->info.dt;
140: dt[1] = ((PC_Factor*)icc)->info.dtcol;
141: dt[2] = ((PC_Factor*)icc)->info.dtcount;
142: PetscInt dtmax = 3;
143: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
144: if (flg) {
145: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
146: }
147: */
148: PetscOptionsTail();
149: return(0);
150: }
152: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
154: /*MC
155: PCICC - Incomplete Cholesky factorization preconditioners.
157: Options Database Keys:
158: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
159: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
160: its factorization (overwrites original matrix)
161: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
162: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
164: Level: beginner
166: Notes:
167: Only implemented for some matrix formats. Not implemented in parallel.
169: For BAIJ matrices this implements a point block ICC.
171: The Manteuffel shift is only implemented for matrices with block size 1
173: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
174: to turn off the shift.
176: References:
177: . 1. - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
178: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
179: Science and Engineering, Kluwer.
182: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
183: PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
184: PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
185: PCFactorSetLevels()
187: M*/
189: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
190: {
192: PC_ICC *icc;
195: PetscNewLog(pc,&icc);
196: pc->data = (void*)icc;
197: PCFactorInitialize(pc);
198: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);
200: ((PC_Factor*)icc)->factortype = MAT_FACTOR_ICC;
201: ((PC_Factor*)icc)->info.fill = 1.0;
202: ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT;
203: ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
205: pc->ops->apply = PCApply_ICC;
206: pc->ops->matapply = PCMatApply_ICC;
207: pc->ops->applytranspose = PCApply_ICC;
208: pc->ops->setup = PCSetUp_ICC;
209: pc->ops->reset = PCReset_ICC;
210: pc->ops->destroy = PCDestroy_ICC;
211: pc->ops->setfromoptions = PCSetFromOptions_ICC;
212: pc->ops->view = PCView_Factor;
213: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
214: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
215: return(0);
216: }