Actual source code: mlInterface.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: mlInterface.c,v 1.1 1999/06/09 20:24:12 knepley Exp $";
3: #endif
4: /*
5: Defines the interface for the multilevel preconditioner due to
6: Vivek Sarin for AIJ matrices.
7: */
8: #include src/sles/pc/pcimpl.h
9: #include petscsnes.h
10: #include ml.h
12: #undef __FUNCT__
14: /*@C PCMultiLevelSetFields
15: This function sets the shape and test function fields for the
16: gradient operator. It must be called prior to PCSetUp().
18: Collective on PC
20: Input Parameters:
21: + pc - The preconditioner context
22: . sField - The shape function field
23: - tField - The test function field
25: Level: intermediate
27: .keywords multilevel, fields
28: .seealso PCMultiLevelSetGradientOperator, PCMultiLevelSetNonlinearIterate
29: @*/
30: int PCMultiLevelSetFields(PC pc, int sField, int tField)
31: {
32: PC_Multilevel *ml;
36: ml = (PC_Multilevel *) pc->data;
37: ml->sField = sField;
38: ml->tField = tField;
39: return(0);
40: }
42: #undef __FUNCT__
44: /*@C PCMultiLevelSetNonlinearIterate
45: This function sets the solution vector used in a Newton iteration. If
46: the inner solve is preconditioned with ML, the current iterate is
47: necessary to perform the reduction of the system.
49: Collective on PC
51: Input Parameters:
52: + pc - The preconditioner context
53: - u - The current iterate in a Newton solve
55: Level: intermediate
57: .keywords multilevel, nonlinear
58: .seealso PCMultiLevelSetFields, PCMultiLevelSetGradientOperator
59: @*/
60: int PCMultiLevelSetNonlinearIterate(PC pc, GVec u)
61: {
62: PC_Multilevel *ml;
67: ml = (PC_Multilevel *) pc->data;
68: ml->nonlinearIterate = u;
69: return(0);
70: }
72: #undef __FUNCT__
74: /*@C PCMultiLevelSetGradientOperator
75: This function sets the operator to use as the gradient and its transpose
76: in the multilevel scheme. It must be called prior to PCSetUp().
78: Collective on PC
80: Input Parameters:
81: + pc - The preconditioner context
82: . grad - The gradient operator
83: . div - The divergence operator
84: - alpha - The scalar multiplier for the gradient operator
86: Level: intermediate
88: .keywords multilevel, operator
89: .seealso PCMultiLevelSetFields, PCMultiLevelSetNonlinearIterate
90: @*/
91: int PCMultiLevelSetGradientOperator(PC pc, int Grad, int Div, PetscScalar alpha)
92: {
93: PC_Multilevel *ml;
97: ml = (PC_Multilevel *) pc->data;
98: ml->gradOp = Grad;
99: ml->divOp = Div;
100: ml->gradAlpha = alpha;
101: return(0);
102: }
104: #undef __FUNCT__
106: /*@C PCMultiLevelBuildSolution
107: This function recovers the solution of the global problem
108: from the solution of the reduced system.
110: Collective on GVec
112: Input Parameters:
113: + pc - The preconditioner context
114: - u - The solution of the reduced system
116: Output Parameter:
117: . u - The global solution
119: Level: beginner
121: .keywords multilevel, solution
122: .seealso PCMultiLevelApplyGradientGetMultiplier
123: @*/
124: int PCMultiLevelBuildSolution(PC pc, GVec u)
125: {
126: PC_Multilevel *ml;
127: Grid grid;
128: PetscScalar alpha;
129: int ierr;
134: if (pc->setupcalled < 2) {
135: PCSetUp(pc);
136: }
138: ml = (PC_Multilevel *) pc->data;
139: if (ml->reduceSol != PETSC_NULL) {
140: ierr = GVecGetGrid(u, &grid);
141: ierr = GridGetBCMultiplier(grid, &alpha);
142: alpha = -alpha;
143: ierr = VecAXPY(&alpha, ml->reduceSol, u);
144: }
145: return(0);
146: }
148: #undef __FUNCT__
150: /*@C PCMultiLevelGetMultiplier
151: This function recovers the multiplier of the global problem
152: from the solution of the reduced system.
154: Collective on GVec
156: Input Parameters:
157: + pc - The preconditioner context
158: - u - The solution of the reduced system
160: Output Parameter:
161: . p - The global multiplier
163: Level: beginner
165: .keywords multilevel, solution, multiplier
166: .seealso PCMultiLevelBuildSolution
167: @*/
168: int PCMultiLevelGetMultiplier(PC pc, GVec u, GVec p)
169: {
170: GVec tempVec;
171: PetscScalar minusOne = -1.0;
172: PetscScalar zero = 0.0;
173: int ierr;
179: if ((pc->vec == PETSC_NULL) || (pc->mat == PETSC_NULL)) {
180: VecSet(&zero, p);
181: return(0);
182: }
183: if (pc->setupcalled < 2) {
184: PCSetUp(pc);
185: }
187: /* Multiply P_2 (P^{-1} u)_2 by A */
188: GVecDuplicate(u, &tempVec);
189: MatMult(pc->mat, u, tempVec);
191: /* Let A P_2 (P^{-1} u)_2 = f - A P_1 D^{-1} Z^T g - A P_2 (P^{-1} u)_2 --
192: This is really unstaisfactory, but how do I get the Rhs after solution of
193: the system. Perhaps I should have a flag for storing the multiplier.
194: */
195: VecAYPX(&minusOne, pc->vec, tempVec);
197: /* Apply P^T_1 */
198: PCMultiLevelApplyP1Trans(pc, tempVec, p);
200: /* Apply Z D^{-1} */
201: PCMultiLevelApplyDInv(pc, p, p);
202: PCMultiLevelApplyV(pc, p, p);
204: /* Cleanup */
205: VecDestroy(tempVec);
206: return(0);
207: }