Actual source code: gsles.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: gsles.c,v 1.5 2000/01/10 03:54:18 knepley Exp $";
3: #endif
5: #include petscsles.h
6: #include gsolver.h
8: #undef __FUNCT__
10: /*@
11: GVecSolutionKSPMonitor - Monitors solution at each KSP iteration.
13: Collective on KSP
15: Input Parameters:
16: . ksp - The Krylov subspace context
17: . it - The number of iterations so far
18: . rnorm - The current (approximate) residual norm
19: . ctx - A viewer
21: Level: advanced
23: .keywords: grid vector, ksp, monitor, solution
24: .seealso: KSPDefaultMonitor(),KSPSetMonitor(),GVecResidualKSPMonitor(),
25: GVecRhsKSPMonitor(), GVecErrorKSPMonitor()
26: @*/
27: int GVecSolutionKSPMonitor(KSP ksp, int it, PetscReal rnorm, void *ctx)
28: {
29: PetscViewer viewer = (PetscViewer) ctx;
30: Vec x;
31: int ierr;
36: KSPBuildSolution(ksp, PETSC_NULL, &x);
37: VecView(x, viewer);
38: return(0);
39: }
41: #undef __FUNCT__
43: /*@
44: GVecResidualKSPMonitor - Monitors residual at each KSP iteration.
46: Collective on KSP
48: Input Parameters:
49: . ksp - The Krylov subspace context
50: . it - The number of iterations so far
51: . rnorm - The current (approximate) residual norm
52: . ctx - A viewer
54: Level: advanced
56: .keywords: grid vector, ksp, monitor, residual
57: .seealso: KSPDefaultMonitor(),KSPSetMonitor(),GVecSolutionKSPMonitor(),
58: GVecRhsKSPMonitor(), GVecErrorKSPMonitor()
59: @*/
60: int GVecResidualKSPMonitor(KSP ksp, int it, PetscReal rnorm, void *ctx)
61: {
62: PetscViewer viewer = (PetscViewer) ctx;
63: Vec r;
64: int ierr;
69: KSPBuildResidual(ksp, PETSC_NULL, PETSC_NULL, &r);
70: VecView(r, viewer);
71: return(0);
72: }
74: #undef __FUNCT__
76: /*@
77: GVecRhsKSPMonitor - Displays the right hand side for the linear system at the first iteration.
79: Collective on KSP
81: Input Parameters:
82: . ksp - The Krylov subspace context
83: . it - The number of iterations so far
84: . rnorm - The current (approximate) residual norm
85: . ctx - A viewer
87: Level: advanced
89: .keywords: grid vector, ksp, monitor, rhs
90: .seealso: KSPDefaultMonitor(),KSPSetMonitor(),GVecSolutionKSPMonitor(),
91: GVecResidualKSPMonitor(), GVecErrorKSPMonitor()
92: @*/
93: int GVecRhsKSPMonitor(KSP ksp, int it, PetscReal rnorm, void *ctx)
94: {
95: PetscViewer viewer = (PetscViewer) ctx;
96: Vec b;
97: int ierr;
102: if (!it) {
103: KSPGetRhs(ksp, &b);
104: VecView(b, viewer);
105: }
106: return(0);
107: }
109: #undef __FUNCT__
111: /*@
112: GVecErrorKSPMonitor - Displays the error at each iteration.
114: Collective on KSP
116: Input Parameters:
117: . ksp - The Krylov subspace context
118: . it - The number of iterations so far
119: . rnorm - The current (approximate) residual norm
120: . errorCtx - A GVecErrorKSPMonitorCtx
122: Notes:
123: The final argument to KSPSetMonitor() with this routine must be
124: a pointer to a GVecErrorKSPMonitorCtx.
126: Level: advanced
128: .keywords: grid vector, ksp, monitor, error
129: .seealso: KSPDefaultMonitor(),KSPSetMonitor(),,GVecSolutionKSPMonitor(),
130: GVecResidualKSPMonitor(), GVecRhsKSPMonitor()
131: @*/
132: int GVecErrorKSPMonitor(KSP ksp, int it, PetscReal rnorm, void *errorCtx)
133: {
134: GVecErrorKSPMonitorCtx *ctx = (GVecErrorKSPMonitorCtx *) errorCtx;
135: PetscScalar mone = -1.0;
136: GVec x, e;
137: PetscReal norm_2, norm_max;
138: int ierr;
142: KSPBuildSolution(ksp, PETSC_NULL, &x);
143: GVecGetWorkGVec(ctx->solution, &e);
144: VecWAXPY(&mone, x, ctx->solution, e);
145: VecView(e, ctx->error_viewer);
147: /* Compute 2-norm and max-norm of error */
148: if (ctx->norm_error_viewer) {
149: VecNorm(e, NORM_2, &norm_2);
150: VecNorm(e, NORM_MAX, &norm_max);
151: if (ctx->isInner == PETSC_FALSE) {
152: PetscViewerASCIIPrintf(ctx->norm_error_viewer, "Iteration %d residual norm %g error 2-norm %g error max-norm %gn",
153: it, rnorm, norm_2, norm_max);
154: } else {
155: PetscViewerASCIIPrintf(ctx->norm_error_viewer, " Inner iteration %d residual norm %g error 2-norm %g error max-norm %gn",
156: it, rnorm, norm_2, norm_max);
157: }
158: }
160: GVecRestoreWorkGVec(ctx->solution, &e);
161: return(0);
162: }
164: #include petscmg.h
165: #undef __FUNCT__
167: /* ------------------------------------------------------------------------------*/
168: /*
169: GVecPCMGSetMonitor - Sets GVec KSP monitors for each level of multigrid.
171: Input Parameters:
172: . pc - preconditioner context for multigrid
173: */
174: int GVecPCMGSetMonitor(PC pc, int views)
175: {
176: char name[25];
177: int levels, i,ierr, solution,rhs,residual,all;
178: SLES subsles;
179: KSP ksp;
180: PetscViewer *viewers;
181: PetscTruth match;
185: PetscTypeCompare((PetscObject) pc, PCMG, &match);
186: if (match == PETSC_FALSE) return(0);
188: solution = (views & GVECPCMGMONITOR_SOLUTION);
189: residual = (views & GVECPCMGMONITOR_RESIDUAL);
190: rhs = (views & GVECPCMGMONITOR_RHS);
191: all = (solution != 0) + (residual != 0) + (rhs != 0);
193: MGGetLevels(pc, &levels);
195: PetscMalloc(all*levels * sizeof(PetscViewer), &viewers);
197: for(i = 0; i < levels; i++) {
198: MGGetSmoother(pc, i, &subsles);
199: SLESGetKSP(subsles, &ksp);
200: if (rhs) {
201: sprintf(name, "Right hand side %d",i);
202: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, name, PETSC_DECIDE, PETSC_DECIDE, 200, 200, viewers);
203: KSPSetMonitor(ksp, GVecRhsKSPMonitor, *viewers++, PETSC_NULL);
204: }
205: if (solution) {
206: sprintf(name, "Solution %d",i);
207: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, name, PETSC_DECIDE, PETSC_DECIDE, 200, 200, viewers);
208: KSPSetMonitor(ksp, GVecSolutionKSPMonitor, *viewers++, PETSC_NULL);
209: }
210: if (residual) {
211: sprintf(name, "Residual %d",i);
212: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, name, PETSC_DECIDE, PETSC_DECIDE, 200, 200, viewers);
213: KSPSetMonitor(ksp, GVecResidualKSPMonitor, *viewers++, PETSC_NULL);
214: }
215: }
217: return(0);
218: }
220: EXTERN_C_BEGIN
221: #undef __FUNCT__
223: int GVecKSPOptionsChecker_Private(KSP ksp)
224: {
225: char *prefix, string[64], *p;
226: int ierr, loc[4], nmax;
227: MPI_Comm comm;
228: PetscViewer viewer;
229: PetscDraw draw;
230: PetscTruth opt;
235: KSPGetOptionsPrefix(ksp, &prefix);
236: PetscObjectGetComm((PetscObject) ksp, &comm);
238: nmax = 4;
239: loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
240: PetscOptionsGetIntArray(prefix, "-gvec_ksp_solutionmonitor", loc, &nmax, &opt);
241: if (opt) {
242: PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);
243: PetscViewerDrawGetDraw(viewer, 0, &draw);
244: PetscDrawSetTitle(draw, "Approx. Solution");
245: PetscLogObjectParent(ksp, (PetscObject) viewer);
246: KSPSetMonitor(ksp, GVecSolutionKSPMonitor, (void *) viewer, PETSC_NULL);
247: }
248: nmax = 4;
249: loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
250: PetscOptionsGetIntArray(prefix, "-gvec_ksp_residualmonitor", loc, &nmax, &opt);
251: if (opt) {
252: PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);
253: PetscViewerDrawGetDraw(viewer, 0, &draw);
254: PetscDrawSetTitle(draw, "Residual");
255: PetscLogObjectParent(ksp, (PetscObject) viewer);
256: KSPSetMonitor(ksp, GVecResidualKSPMonitor, (void *) viewer, PETSC_NULL);
257: }
258: nmax = 4;
259: loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
260: PetscOptionsGetIntArray(prefix, "-gvec_ksp_rhsmonitor", loc, &nmax, &opt);
261: if (opt) {
262: PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);
263: PetscViewerDrawGetDraw(viewer, 0, &draw);
264: PetscDrawSetTitle(draw, "Rhs");
265: PetscLogObjectParent(ksp, (PetscObject) viewer);
266: KSPSetMonitor(ksp, GVecRhsKSPMonitor, (void *) viewer, PETSC_NULL);
267: }
268: PetscOptionsGetString(prefix, "-gvec_mg_monitor", string, 64, &opt);
269: if (opt) {
270: PC pc;
271: int all = 0;
273: PetscStrstr(string, "nosolution", &p);
274: if (!p) all |= GVECPCMGMONITOR_SOLUTION;
275: PetscStrstr(string, "noresidual", &p);
276: if (!p) all |= GVECPCMGMONITOR_RESIDUAL;
277: PetscStrstr(string, "norhs", &p);
278: if (!p) all |= GVECPCMGMONITOR_RHS;
279: KSPGetPC(ksp, &pc);
280: GVecPCMGSetMonitor(pc, all);
281: }
282: PetscOptionsHasName(PETSC_NULL, "-help", &opt);
283: if (opt) {
284: char *pprefix;
285: int len = 2;
287: if (prefix != PETSC_NULL) {
288: PetscStrlen(prefix, &len);
289: }
290: PetscMalloc((len+2) * sizeof(char), &pprefix);
291: PetscStrcpy(pprefix, "-");
292: if (prefix != PETSC_NULL)
293: PetscStrcat(pprefix, prefix);
294: PetscPrintf(comm," Additional KSP Monitor options for grid vectorsn");
295: PetscPrintf(comm," %sgvec_ksp_solutionmonitorn", pprefix);
296: PetscPrintf(comm," %sgvec_ksp_residualmonitorn", pprefix);
297: PetscPrintf(comm," %sgvec_ksp_rhsmonitorn", pprefix);
298: PetscPrintf(comm," %sgvec_mg_monitor [nosolution,norhs,noresidual]n", pprefix);
299: PetscFree(pprefix);
300: }
301: return(0);
302: }
303: EXTERN_C_END
305: /*--------------------------------------------------- Uzawa Methods --------------------------------------------------*/
306: #undef __FUNCT__
308: /*@
309: GMatCreateUzawa - Creates a container matrix for the Uzawa system matrix B^T A^{-1} B.
311: Input Parameter:
312: + A - The square (0,0) block of the original saddle-point matrix
313: . B - The rectangular (0,1) block of the original saddle-point matrix
314: - sles - The solver used to invert A
316: Output Parameter:
317: . gmat - The discrete operator
319: Level: advanced
321: .seealso: GVecCreate()
322: @*/
323: int GMatCreateUzawa(GMat A, GMat B, SLES sles, GMat *gmat)
324: {
325: UzawaContext *uCtx;
326: MPI_Comm comm;
327: Grid grid, grid2;
328: int M, m, N, n;
329: int O, o, P, p;
330: int ierr;
337: GMatGetGrid(A, &grid);
338: GMatGetGrid(B, &grid2);
341: if (grid != grid2) SETERRQ(PETSC_ERR_ARG_INCOMP, "Matrices must all have the same underlying grid.");
342: PetscNew(UzawaContext, &uCtx);
343: uCtx->A = A;
344: uCtx->B = B;
345: uCtx->sles = sles;
346: MatGetSize(A, &M, &N);
347: MatGetLocalSize(A, &m, &n);
348: MatGetSize(B, &O, &P);
349: MatGetLocalSize(B, &o, &p);
350: if ((m != o) || (n != o) || (M != O) || (N != O)) {
351: SETERRQ(PETSC_ERR_ARG_INCOMP, "Incommensurate sizes, B^T A B must be a legal matrix");
352: }
353: PetscObjectGetComm((PetscObject) grid, &comm);
354: MatCreateShell(comm, p, p, P, P, uCtx, gmat);
355: MatShellSetOperation(*gmat, MATOP_MULT, (void (*)(void)) GMatMatMultUzawa);
356: VecCreateMPI(comm, m, M, &uCtx->work);
357: VecDuplicate(uCtx->work, &uCtx->work2);
358: return(0);
359: }
361: #undef __FUNCT__
363: /*@
364: GMatMatMultUzawa - This function applies B^T A^{-1} B to a vector,
365: which is the Schur complement matrix.
367: Input Parameters:
368: + mat - The grid matrix
369: - x - The input grid vector
371: Output Parameter:
372: . y - The output grid vector B^T A^{-1} B x
374: Level: advanced
376: .seealso: GMatEvaluateOperatorGalerkin
377: @*/
378: int GMatMatMultUzawa(GMat mat, GVec x, GVec y)
379: {
380: UzawaContext *ctx;
381: int its;
382: int ierr;
388: MatShellGetContext(mat, (void **) &ctx);
389: MatMult(ctx->B, x, ctx->work);
390: SLESSolve(ctx->sles, ctx->work, ctx->work2, &its);
391: MatMultTranspose(ctx->B, ctx->work2, y);
392: return(0);
393: }
395: #undef __FUNCT__
397: /*@
398: GMatDestroyUzawa - Destroys a container matrix for the Uzawa
399: system matrix B^T A B.
401: Input Parameter:
402: . gmat - The container matrix from GMatCreateUzawa()
404: Level: advanced
406: .seealso: GVecCreate()
407: @*/
408: int GMatDestroyUzawa(GMat gmat)
409: {
410: UzawaContext *ctx;
411: int ierr;
416: MatShellGetContext(gmat, (void **) &ctx);
417: if (ctx != PETSC_NULL) {
418: VecDestroy(ctx->work);
419: VecDestroy(ctx->work2);
420: PetscFree(ctx);
421: }
422: MatDestroy(gmat);
423: return(0);
424: }