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