Actual source code: ex3.c
1: static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
2: This example employs a user-defined monitoring routine and optionally a user-defined\n\
3: routine to check candidate iterates produced by line search routines.\n\
4: The command line options include:\n\
5: -pre_check_iterates : activate checking of iterates\n\
6: -post_check_iterates : activate checking of iterates\n\
7: -check_tol <tol>: set tolerance for iterate checking\n\
8: -user_precond : activate a (trivial) user-defined preconditioner\n\n";
10: /*T
11: Concepts: SNES^basic parallel example
12: Concepts: SNES^setting a user-defined monitoring routine
13: Processors: n
14: T*/
18: /*
19: Include "petscdm.h" so that we can use data management objects (DMs)
20: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
21: Include "petscsnes.h" so that we can use SNES solvers. Note that this
22: file automatically includes:
23: petscsys.h - base PETSc routines
24: petscvec.h - vectors
25: petscmat.h - matrices
26: petscis.h - index sets
27: petscksp.h - Krylov subspace methods
28: petscviewer.h - viewers
29: petscpc.h - preconditioners
30: petscksp.h - linear solvers
31: */
33: #include <petscdm.h>
34: #include <petscdmda.h>
35: #include <petscsnes.h>
37: /*
38: User-defined routines.
39: */
40: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
41: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
42: PetscErrorCode FormInitialGuess(Vec);
43: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
44: PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
45: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
46: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
47: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
49: /*
50: User-defined application context
51: */
52: typedef struct {
53: DM da; /* distributed array */
54: Vec F; /* right-hand-side of PDE */
55: PetscMPIInt rank; /* rank of processor */
56: PetscMPIInt size; /* size of communicator */
57: PetscReal h; /* mesh spacing */
58: PetscBool sjerr; /* if or not to test jacobian domain error */
59: } ApplicationCtx;
61: /*
62: User-defined context for monitoring
63: */
64: typedef struct {
65: PetscViewer viewer;
66: } MonitorCtx;
68: /*
69: User-defined context for checking candidate iterates that are
70: determined by line search methods
71: */
72: typedef struct {
73: Vec last_step; /* previous iterate */
74: PetscReal tolerance; /* tolerance for changes between successive iterates */
75: ApplicationCtx *user;
76: } StepCheckCtx;
78: typedef struct {
79: PetscInt its0; /* num of prevous outer KSP iterations */
80: } SetSubKSPCtx;
82: int main(int argc,char **argv)
83: {
84: SNES snes; /* SNES context */
85: SNESLineSearch linesearch; /* SNESLineSearch context */
86: Mat J; /* Jacobian matrix */
87: ApplicationCtx ctx; /* user-defined context */
88: Vec x,r,U,F; /* vectors */
89: MonitorCtx monP; /* monitoring context */
90: StepCheckCtx checkP; /* step-checking context */
91: SetSubKSPCtx checkP1;
92: PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
93: PetscScalar xp,*FF,*UU,none = -1.0;
95: PetscInt its,N = 5,i,maxit,maxf,xs,xm;
96: PetscReal abstol,rtol,stol,norm;
97: PetscBool flg,viewinitial = PETSC_FALSE;
100: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
101: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
102: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
103: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
104: ctx.h = 1.0/(N-1);
105: ctx.sjerr = PETSC_FALSE;
106: PetscOptionsGetBool(NULL,NULL,"-test_jacobian_domain_error",&ctx.sjerr,NULL);
107: PetscOptionsGetBool(NULL,NULL,"-view_initial",&viewinitial,NULL);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Create nonlinear solver context
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: SNESCreate(PETSC_COMM_WORLD,&snes);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Create vector data structures; set function evaluation routine
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: /*
120: Create distributed array (DMDA) to manage parallel grid and vectors
121: */
122: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);
123: DMSetFromOptions(ctx.da);
124: DMSetUp(ctx.da);
126: /*
127: Extract global and local vectors from DMDA; then duplicate for remaining
128: vectors that are the same types
129: */
130: DMCreateGlobalVector(ctx.da,&x);
131: VecDuplicate(x,&r);
132: VecDuplicate(x,&F); ctx.F = F;
133: VecDuplicate(x,&U);
135: /*
136: Set function evaluation routine and vector. Whenever the nonlinear
137: solver needs to compute the nonlinear function, it will call this
138: routine.
139: - Note that the final routine argument is the user-defined
140: context that provides application-specific data for the
141: function evaluation routine.
142: */
143: SNESSetFunction(snes,r,FormFunction,&ctx);
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Create matrix data structure; set Jacobian evaluation routine
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: MatCreate(PETSC_COMM_WORLD,&J);
150: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
151: MatSetFromOptions(J);
152: MatSeqAIJSetPreallocation(J,3,NULL);
153: MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);
155: /*
156: Set Jacobian matrix data structure and default Jacobian evaluation
157: routine. Whenever the nonlinear solver needs to compute the
158: Jacobian matrix, it will call this routine.
159: - Note that the final routine argument is the user-defined
160: context that provides application-specific data for the
161: Jacobian evaluation routine.
162: */
163: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
165: /*
166: Optionally allow user-provided preconditioner
167: */
168: PetscOptionsHasName(NULL,NULL,"-user_precond",&flg);
169: if (flg) {
170: KSP ksp;
171: PC pc;
172: SNESGetKSP(snes,&ksp);
173: KSPGetPC(ksp,&pc);
174: PCSetType(pc,PCSHELL);
175: PCShellSetApply(pc,MatrixFreePreconditioner);
176: }
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Customize nonlinear solver; set runtime options
180: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182: /*
183: Set an optional user-defined monitoring routine
184: */
185: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
186: SNESMonitorSet(snes,Monitor,&monP,0);
188: /*
189: Set names for some vectors to facilitate monitoring (optional)
190: */
191: PetscObjectSetName((PetscObject)x,"Approximate Solution");
192: PetscObjectSetName((PetscObject)U,"Exact Solution");
194: /*
195: Set SNES/KSP/KSP/PC runtime options, e.g.,
196: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
197: */
198: SNESSetFromOptions(snes);
200: /*
201: Set an optional user-defined routine to check the validity of candidate
202: iterates that are determined by line search methods
203: */
204: SNESGetLineSearch(snes, &linesearch);
205: PetscOptionsHasName(NULL,NULL,"-post_check_iterates",&post_check);
207: if (post_check) {
208: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
209: SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
210: VecDuplicate(x,&(checkP.last_step));
212: checkP.tolerance = 1.0;
213: checkP.user = &ctx;
215: PetscOptionsGetReal(NULL,NULL,"-check_tol",&checkP.tolerance,NULL);
216: }
218: PetscOptionsHasName(NULL,NULL,"-post_setsubksp",&post_setsubksp);
219: if (post_setsubksp) {
220: PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
221: SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
222: }
224: PetscOptionsHasName(NULL,NULL,"-pre_check_iterates",&pre_check);
225: if (pre_check) {
226: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
227: SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
228: }
230: /*
231: Print parameters used for convergence testing (optional) ... just
232: to demonstrate this routine; this information is also printed with
233: the option -snes_view
234: */
235: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
236: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);
238: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: Initialize application:
240: Store right-hand-side of PDE and exact solution
241: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
243: /*
244: Get local grid boundaries (for 1-dimensional DMDA):
245: xs, xm - starting grid index, width of local grid (no ghost points)
246: */
247: DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
249: /*
250: Get pointers to vector data
251: */
252: DMDAVecGetArray(ctx.da,F,&FF);
253: DMDAVecGetArray(ctx.da,U,&UU);
255: /*
256: Compute local vector entries
257: */
258: xp = ctx.h*xs;
259: for (i=xs; i<xs+xm; i++) {
260: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
261: UU[i] = xp*xp*xp;
262: xp += ctx.h;
263: }
265: /*
266: Restore vectors
267: */
268: DMDAVecRestoreArray(ctx.da,F,&FF);
269: DMDAVecRestoreArray(ctx.da,U,&UU);
270: if (viewinitial) {
271: VecView(U,0);
272: VecView(F,0);
273: }
275: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276: Evaluate initial guess; then solve nonlinear system
277: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279: /*
280: Note: The user should initialize the vector, x, with the initial guess
281: for the nonlinear solver prior to calling SNESSolve(). In particular,
282: to employ an initial guess of zero, the user should explicitly set
283: this vector to zero by calling VecSet().
284: */
285: FormInitialGuess(x);
286: SNESSolve(snes,NULL,x);
287: SNESGetIterationNumber(snes,&its);
288: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
290: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291: Check solution and clean up
292: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
294: /*
295: Check the error
296: */
297: VecAXPY(x,none,U);
298: VecNorm(x,NORM_2,&norm);
299: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its);
300: if (ctx.sjerr) {
301: SNESType snestype;
302: SNESGetType(snes,&snestype);
303: PetscPrintf(PETSC_COMM_WORLD,"SNES Type %s\n",snestype);
304: }
306: /*
307: Free work space. All PETSc objects should be destroyed when they
308: are no longer needed.
309: */
310: PetscViewerDestroy(&monP.viewer);
311: if (post_check) {VecDestroy(&checkP.last_step);}
312: VecDestroy(&x);
313: VecDestroy(&r);
314: VecDestroy(&U);
315: VecDestroy(&F);
316: MatDestroy(&J);
317: SNESDestroy(&snes);
318: DMDestroy(&ctx.da);
319: PetscFinalize();
320: return ierr;
321: }
323: /* ------------------------------------------------------------------- */
324: /*
325: FormInitialGuess - Computes initial guess.
327: Input/Output Parameter:
328: . x - the solution vector
329: */
330: PetscErrorCode FormInitialGuess(Vec x)
331: {
333: PetscScalar pfive = .50;
336: VecSet(x,pfive);
337: return(0);
338: }
340: /* ------------------------------------------------------------------- */
341: /*
342: FormFunction - Evaluates nonlinear function, F(x).
344: Input Parameters:
345: . snes - the SNES context
346: . x - input vector
347: . ctx - optional user-defined context, as set by SNESSetFunction()
349: Output Parameter:
350: . f - function vector
352: Note:
353: The user-defined context can contain any application-specific
354: data needed for the function evaluation.
355: */
356: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
357: {
358: ApplicationCtx *user = (ApplicationCtx*) ctx;
359: DM da = user->da;
360: PetscScalar *xx,*ff,*FF,d;
362: PetscInt i,M,xs,xm;
363: Vec xlocal;
366: DMGetLocalVector(da,&xlocal);
367: /*
368: Scatter ghost points to local vector, using the 2-step process
369: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
370: By placing code between these two statements, computations can
371: be done while messages are in transition.
372: */
373: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
374: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
376: /*
377: Get pointers to vector data.
378: - The vector xlocal includes ghost point; the vectors x and f do
379: NOT include ghost points.
380: - Using DMDAVecGetArray() allows accessing the values using global ordering
381: */
382: DMDAVecGetArray(da,xlocal,&xx);
383: DMDAVecGetArray(da,f,&ff);
384: DMDAVecGetArray(da,user->F,&FF);
386: /*
387: Get local grid boundaries (for 1-dimensional DMDA):
388: xs, xm - starting grid index, width of local grid (no ghost points)
389: */
390: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
391: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
393: /*
394: Set function values for boundary points; define local interior grid point range:
395: xsi - starting interior grid index
396: xei - ending interior grid index
397: */
398: if (xs == 0) { /* left boundary */
399: ff[0] = xx[0];
400: xs++;xm--;
401: }
402: if (xs+xm == M) { /* right boundary */
403: ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
404: xm--;
405: }
407: /*
408: Compute function over locally owned part of the grid (interior points only)
409: */
410: d = 1.0/(user->h*user->h);
411: for (i=xs; i<xs+xm; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
413: /*
414: Restore vectors
415: */
416: DMDAVecRestoreArray(da,xlocal,&xx);
417: DMDAVecRestoreArray(da,f,&ff);
418: DMDAVecRestoreArray(da,user->F,&FF);
419: DMRestoreLocalVector(da,&xlocal);
420: return(0);
421: }
423: /* ------------------------------------------------------------------- */
424: /*
425: FormJacobian - Evaluates Jacobian matrix.
427: Input Parameters:
428: . snes - the SNES context
429: . x - input vector
430: . dummy - optional user-defined context (not used here)
432: Output Parameters:
433: . jac - Jacobian matrix
434: . B - optionally different preconditioning matrix
435: . flag - flag indicating matrix structure
436: */
437: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
438: {
439: ApplicationCtx *user = (ApplicationCtx*) ctx;
440: PetscScalar *xx,d,A[3];
442: PetscInt i,j[3],M,xs,xm;
443: DM da = user->da;
446: if (user->sjerr) {
447: SNESSetJacobianDomainError(snes);
448: return(0);
449: }
450: /*
451: Get pointer to vector data
452: */
453: DMDAVecGetArrayRead(da,x,&xx);
454: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
456: /*
457: Get range of locally owned matrix
458: */
459: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
461: /*
462: Determine starting and ending local indices for interior grid points.
463: Set Jacobian entries for boundary points.
464: */
466: if (xs == 0) { /* left boundary */
467: i = 0; A[0] = 1.0;
469: MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
470: xs++;xm--;
471: }
472: if (xs+xm == M) { /* right boundary */
473: i = M-1;
474: A[0] = 1.0;
475: MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
476: xm--;
477: }
479: /*
480: Interior grid points
481: - Note that in this case we set all elements for a particular
482: row at once.
483: */
484: d = 1.0/(user->h*user->h);
485: for (i=xs; i<xs+xm; i++) {
486: j[0] = i - 1; j[1] = i; j[2] = i + 1;
487: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
488: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
489: }
491: /*
492: Assemble matrix, using the 2-step process:
493: MatAssemblyBegin(), MatAssemblyEnd().
494: By placing code between these two statements, computations can be
495: done while messages are in transition.
497: Also, restore vector.
498: */
500: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
501: DMDAVecRestoreArrayRead(da,x,&xx);
502: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
504: return(0);
505: }
507: /* ------------------------------------------------------------------- */
508: /*
509: Monitor - Optional user-defined monitoring routine that views the
510: current iterate with an x-window plot. Set by SNESMonitorSet().
512: Input Parameters:
513: snes - the SNES context
514: its - iteration number
515: norm - 2-norm function value (may be estimated)
516: ctx - optional user-defined context for private data for the
517: monitor routine, as set by SNESMonitorSet()
519: Note:
520: See the manpage for PetscViewerDrawOpen() for useful runtime options,
521: such as -nox to deactivate all x-window output.
522: */
523: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
524: {
526: MonitorCtx *monP = (MonitorCtx*) ctx;
527: Vec x;
530: PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);
531: SNESGetSolution(snes,&x);
532: VecView(x,monP->viewer);
533: return(0);
534: }
536: /* ------------------------------------------------------------------- */
537: /*
538: PreCheck - Optional user-defined routine that checks the validity of
539: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
541: Input Parameters:
542: snes - the SNES context
543: xcurrent - current solution
544: y - search direction and length
546: Output Parameters:
547: y - proposed step (search direction and length) (possibly changed)
548: changed_y - tells if the step has changed or not
549: */
550: PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
551: {
553: *changed_y = PETSC_FALSE;
554: return(0);
555: }
557: /* ------------------------------------------------------------------- */
558: /*
559: PostCheck - Optional user-defined routine that checks the validity of
560: candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
562: Input Parameters:
563: snes - the SNES context
564: ctx - optional user-defined context for private data for the
565: monitor routine, as set by SNESLineSearchSetPostCheck()
566: xcurrent - current solution
567: y - search direction and length
568: x - the new candidate iterate
570: Output Parameters:
571: y - proposed step (search direction and length) (possibly changed)
572: x - current iterate (possibly modified)
574: */
575: PetscErrorCode PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
576: {
578: PetscInt i,iter,xs,xm;
579: StepCheckCtx *check;
580: ApplicationCtx *user;
581: PetscScalar *xa,*xa_last,tmp;
582: PetscReal rdiff;
583: DM da;
584: SNES snes;
587: *changed_x = PETSC_FALSE;
588: *changed_y = PETSC_FALSE;
590: SNESLineSearchGetSNES(linesearch, &snes);
591: check = (StepCheckCtx*)ctx;
592: user = check->user;
593: SNESGetIterationNumber(snes,&iter);
595: /* iteration 1 indicates we are working on the second iteration */
596: if (iter > 0) {
597: da = user->da;
598: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %g\n",iter,(double)check->tolerance);
600: /* Access local array data */
601: DMDAVecGetArray(da,check->last_step,&xa_last);
602: DMDAVecGetArray(da,x,&xa);
603: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
605: /*
606: If we fail the user-defined check for validity of the candidate iterate,
607: then modify the iterate as we like. (Note that the particular modification
608: below is intended simply to demonstrate how to manipulate this data, not
609: as a meaningful or appropriate choice.)
610: */
611: for (i=xs; i<xs+xm; i++) {
612: if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance;
613: else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
614: if (rdiff > check->tolerance) {
615: tmp = xa[i];
616: xa[i] = .5*(xa[i] + xa_last[i]);
617: *changed_x = PETSC_TRUE;
618: PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%g, x_last=%g, diff=%g, x_new=%g\n",i,(double)PetscAbsScalar(tmp),(double)PetscAbsScalar(xa_last[i]),(double)rdiff,(double)PetscAbsScalar(xa[i]));
619: }
620: }
621: DMDAVecRestoreArray(da,check->last_step,&xa_last);
622: DMDAVecRestoreArray(da,x,&xa);
623: }
624: VecCopy(x,check->last_step);
625: return(0);
626: }
628: /* ------------------------------------------------------------------- */
629: /*
630: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
631: e.g,
632: mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
633: Set by SNESLineSearchSetPostCheck().
635: Input Parameters:
636: linesearch - the LineSearch context
637: xcurrent - current solution
638: y - search direction and length
639: x - the new candidate iterate
641: Output Parameters:
642: y - proposed step (search direction and length) (possibly changed)
643: x - current iterate (possibly modified)
645: */
646: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
647: {
649: SetSubKSPCtx *check;
650: PetscInt iter,its,sub_its,maxit;
651: KSP ksp,sub_ksp,*sub_ksps;
652: PC pc;
653: PetscReal ksp_ratio;
654: SNES snes;
657: SNESLineSearchGetSNES(linesearch, &snes);
658: check = (SetSubKSPCtx*)ctx;
659: SNESGetIterationNumber(snes,&iter);
660: SNESGetKSP(snes,&ksp);
661: KSPGetPC(ksp,&pc);
662: PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);
663: sub_ksp = sub_ksps[0];
664: KSPGetIterationNumber(ksp,&its); /* outer KSP iteration number */
665: KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */
667: if (iter) {
668: PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %D %D, subksp_it %D\n",iter,check->its0,its,sub_its);
669: ksp_ratio = ((PetscReal)(its))/check->its0;
670: maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
671: if (maxit < 2) maxit = 2;
672: KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
673: PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %D\n\n",(double)ksp_ratio,maxit);
674: }
675: check->its0 = its; /* save current outer KSP iteration number */
676: return(0);
677: }
679: /* ------------------------------------------------------------------- */
680: /*
681: MatrixFreePreconditioner - This routine demonstrates the use of a
682: user-provided preconditioner. This code implements just the null
683: preconditioner, which of course is not recommended for general use.
685: Input Parameters:
686: + pc - preconditioner
687: - x - input vector
689: Output Parameter:
690: . y - preconditioned vector
691: */
692: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
693: {
695: VecCopy(x,y);
696: return 0;
697: }
700: /*TEST
702: test:
703: args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
705: test:
706: suffix: 2
707: nsize: 3
708: args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
710: test:
711: suffix: 3
712: nsize: 2
713: args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
715: test:
716: suffix: 4
717: args: -nox -pre_check_iterates -post_check_iterates
719: test:
720: suffix: 5
721: requires: double !complex !single
722: nsize: 2
723: args: -nox -snes_test_jacobian -snes_test_jacobian_view
725: test:
726: suffix: 6
727: requires: double !complex !single
728: nsize: 4
729: args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1
731: test:
732: suffix: 7
733: requires: double !complex !single
734: nsize: 4
735: args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontr -snes_check_jacobian_domain_error 1
737: test:
738: suffix: 8
739: requires: double !complex !single
740: nsize: 4
741: args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonrsls -snes_check_jacobian_domain_error 1
743: test:
744: suffix: 9
745: requires: double !complex !single
746: nsize: 4
747: args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonssls -snes_check_jacobian_domain_error 1
749: test:
750: suffix: 10
751: requires: double !complex !single
752: nsize: 4
753: args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1
755: test:
756: suffix: 11
757: requires: double !complex !single
758: nsize: 4
759: args: -test_jacobian_domain_error -snes_converged_reason -snes_type ms -snes_ms_type m62 -snes_ms_damping 0.9 -snes_check_jacobian_domain_error 1
761: test:
762: suffix: 12
763: args: -view_initial
764: filter: grep -v "type:"
766: TEST*/