Actual source code: ex38.c
petsc-3.3-p7 2013-05-11
2: static char help[] = "u`` + u^{2} = f. Different matrices for the Jacobian and the preconditioner.\n\
3: Demonstrates the use of matrix-free Newton-Krylov methods in conjunction\n\
4: with a user-provided preconditioner. Input arguments are:\n\
5: -snes_mf : Use matrix-free Newton methods\n\
6: -user_precond : Employ a user-defined preconditioner. Used only with\n\
7: matrix-free methods in this example.\n\n";
8: /*
9: Modified from ex6.c by Mike McCourt <mccomic@iit.edu>
10: for testing SNESLineSearchSet()
11: */
13: /*T
14: Concepts: SNES^different matrices for the Jacobian and preconditioner;
15: Concepts: SNES^matrix-free methods
16: Concepts: SNES^user-provided preconditioner;
17: Concepts: matrix-free methods
18: Concepts: user-provided preconditioner;
19: Processors: 1
20: T*/
22: /*
23: Include "petscsnes.h" so that we can use SNES solvers. Note that this
24: file automatically includes:
25: petscsys.h - base PETSc routines petscvec.h - vectors
26: petscmat.h - matrices
27: petscis.h - index sets petscksp.h - Krylov subspace methods
28: petscviewer.h - viewers petscpc.h - preconditioners
29: petscksp.h - linear solvers
30: */
31: #include <iostream>
32: using namespace std;
33: #include <petscsnes.h>
35: struct AppCtx{int testint;};
37: /*
38: User-defined routines
39: */
40: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
41: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
42: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
43: PetscErrorCode FormLineSearch(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscBool *);
45: int main(int argc,char **argv)
46: {
47: SNES snes; /* SNES context */
48: KSP ksp; /* KSP context */
49: PC pc; /* PC context */
50: Vec x,r,F; /* vectors */
51: Mat J,JPrec; /* Jacobian,preconditioner matrices */
53: PetscInt it,n = 5,i;
54: PetscMPIInt size;
55: PetscInt *Shistit = 0,Khistl = 200,Shistl = 10;
56: PetscReal h,xp = 0.0,*Khist = 0,*Shist = 0;
57: PetscScalar v,pfive = .5;
58: PetscBool flg;
59: AppCtx user;
61: PetscInitialize(&argc,&argv,(char *)0,help);
62: MPI_Comm_size(PETSC_COMM_WORLD,&size);
63: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
64: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
65: h = 1.0/(n-1);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create nonlinear solver context
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: SNESCreate(PETSC_COMM_WORLD,&snes);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create vector data structures; set function evaluation routine
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: VecCreate(PETSC_COMM_SELF,&x);
78: VecSetSizes(x,PETSC_DECIDE,n);
79: VecSetFromOptions(x);
80: VecDuplicate(x,&r);
81: VecDuplicate(x,&F);
83: SNESSetFunction(snes,r,FormFunction,(void*)F);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Create matrix data structures; set Jacobian evaluation routine
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&J);
90: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,1,PETSC_NULL,&JPrec);
92: /*
93: Note that in this case we create separate matrices for the Jacobian
94: and preconditioner matrix. Both of these are computed in the
95: routine FormJacobian()
96: */
97: SNESSetJacobian(snes,J,JPrec,FormJacobian,0);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Customize nonlinear solver; set runtime options
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: /* Set preconditioner for matrix-free method */
104: flg = PETSC_FALSE;
105: PetscOptionsGetBool(PETSC_NULL,"-snes_mf",&flg,PETSC_NULL);
106: if (flg) {
107: SNESGetKSP(snes,&ksp);
108: KSPGetPC(ksp,&pc);
109: PetscOptionsHasName(PETSC_NULL,"-user_precond",&flg);
110: if (flg) { /* user-defined precond */
111: PCSetType(pc,PCSHELL);
112: PCShellSetApply(pc,MatrixFreePreconditioner);
113: } else {PCSetType(pc,PCNONE);}
114: }
116: SNESSetFromOptions(snes);
117: user.testint = 0;
118: SNESLineSearchSet(snes,FormLineSearch,&user);
120: /*
121: Save all the linear residuals for all the Newton steps; this enables us
122: to retain complete convergence history for printing after the conclusion
123: of SNESSolve(). Alternatively, one could use the monitoring options
124: -snes_monitor -ksp_monitor
125: to see this information during the solver's execution; however, such
126: output during the run distorts performance evaluation data. So, the
127: following is a good option when monitoring code performance, for example
128: when using -log_summary.
129: */
130: PetscOptionsHasName(PETSC_NULL,"-rhistory",&flg);
131: if (flg) {
132: SNESGetKSP(snes,&ksp);
133: PetscMalloc(Khistl*sizeof(PetscReal),&Khist);
134: KSPSetResidualHistory(ksp,Khist,Khistl,PETSC_FALSE);
135: PetscMalloc(Shistl*sizeof(PetscReal),&Shist);
136: PetscMalloc(Shistl*sizeof(PetscInt),&Shistit);
137: SNESSetConvergenceHistory(snes,Shist,Shistit,Shistl,PETSC_FALSE);
138: }
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Initialize application:
142: Store right-hand-side of PDE and exact solution
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: xp = 0.0;
146: for (i=0; i<n; i++) {
147: v = 6.0*xp + pow(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
148: VecSetValues(F,1,&i,&v,INSERT_VALUES);
149: xp += h;
150: }
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Evaluate initial guess; then solve nonlinear system
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: VecSet(x,pfive);
157: SNESSolve(snes,PETSC_NULL,x);
158: SNESGetIterationNumber(snes,&it);
159: PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it);
161: PetscOptionsHasName(PETSC_NULL,"-rhistory",&flg);
162: if (flg) {
163: KSPGetResidualHistory(ksp,PETSC_NULL,&Khistl);
164: PetscRealView(Khistl,Khist,PETSC_VIEWER_STDOUT_SELF);
165: PetscFree(Khist);
166: SNESGetConvergenceHistory(snes,PETSC_NULL,PETSC_NULL,&Shistl);
167: PetscRealView(Shistl,Shist,PETSC_VIEWER_STDOUT_SELF);
168: PetscIntView(Shistl,Shistit,PETSC_VIEWER_STDOUT_SELF);
169: PetscFree(Shist);
170: PetscFree(Shistit);
171: }
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Free work space. All PETSc objects should be destroyed when they
175: are no longer needed.
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178: VecDestroy(&x); VecDestroy(&r);
179: VecDestroy(&F); MatDestroy(&J);
180: MatDestroy(&JPrec); SNESDestroy(&snes);
181: PetscFinalize();
183: return 0;
184: }
186: PetscErrorCode FormLineSearch(SNES snes,void* user,Vec X,Vec F,Vec G,Vec Y,Vec W,PetscReal fnorm,
187: PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag)
188: {
190: PetscScalar mone=-1.0;
191: AppCtx *myguy = (AppCtx*)user;
192: *flag=PETSC_TRUE;
194: cout << "Inside FormLineSearch \n user.testint=" << myguy->testint << endl;
195: ierr=VecNorm(Y,NORM_2,ynorm);
196: ierr=VecWAXPY(W,mone,Y,X); /* W = -Y + X */
197: ierr=SNESComputeFunction(snes,W,G);
198: ierr=VecNorm(G,NORM_2,gnorm);
199: myguy->testint++;
200: return ierr;
201: }
203: /* ------------------------------------------------------------------- */
204: /*
205: FormInitialGuess - Forms initial approximation.
207: Input Parameters:
208: user - user-defined application context
209: X - vector
211: Output Parameter:
212: X - vector
213: */
214: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
215: {
216: PetscScalar *xx,*ff,*FF,d;
218: PetscInt i,n;
220: VecGetArray(x,&xx);
221: VecGetArray(f,&ff);
222: VecGetArray((Vec)dummy,&FF);
223: VecGetSize(x,&n);
224: d = (PetscReal)(n - 1); d = d*d;
225: ff[0] = xx[0];
226: for (i=1; i<n-1; i++) {
227: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
228: }
229: ff[n-1] = xx[n-1] - 1.0;
230: VecRestoreArray(x,&xx);
231: VecRestoreArray(f,&ff);
232: VecRestoreArray((Vec)dummy,&FF);
233: return 0;
234: }
235: /* ------------------------------------------------------------------- */
236: /*
237: FormJacobian - This routine demonstrates the use of different
238: matrices for the Jacobian and preconditioner
240: Input Parameters:
241: . snes - the SNES context
242: . x - input vector
243: . ptr - optional user-defined context, as set by SNESSetJacobian()
245: Output Parameters:
246: . A - Jacobian matrix
247: . B - different preconditioning matrix
248: . flag - flag indicating matrix structure
249: */
250: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *prejac,MatStructure *flag,void *dummy)
251: {
252: PetscScalar *xx,A[3],d;
253: PetscInt i,n,j[3];
256: VecGetArray(x,&xx);
257: VecGetSize(x,&n);
258: d = (PetscReal)(n - 1); d = d*d;
260: /* Form Jacobian. Also form a different preconditioning matrix that
261: has only the diagonal elements. */
262: i = 0; A[0] = 1.0;
263: MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
264: MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
265: for (i=1; i<n-1; i++) {
266: j[0] = i - 1; j[1] = i; j[2] = i + 1;
267: A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d;
268: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
269: MatSetValues(*prejac,1,&i,1,&i,&A[1],INSERT_VALUES);
270: }
271: i = n-1; A[0] = 1.0;
272: MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
273: MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
275: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
276: MatAssemblyBegin(*prejac,MAT_FINAL_ASSEMBLY);
277: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
278: MatAssemblyEnd(*prejac,MAT_FINAL_ASSEMBLY);
280: VecRestoreArray(x,&xx);
281: *flag = SAME_NONZERO_PATTERN;
282: return 0;
283: }
284: /* ------------------------------------------------------------------- */
285: /*
286: MatrixFreePreconditioner - This routine demonstrates the use of a
287: user-provided preconditioner. This code implements just the null
288: preconditioner, which of course is not recommended for general use.
290: Input Parameters:
291: + pc - preconditioner object
292: - x - input vector
294: Output Parameter:
295: . y - preconditioned vector
296: */
297: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
298: {
300: VecCopy(x,y);
301: return 0;
302: }