Actual source code: ex6.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";
9: /*T
10: Concepts: SNES^different matrices for the Jacobian and preconditioner;
11: Concepts: SNES^matrix-free methods
12: Concepts: SNES^user-provided preconditioner;
13: Concepts: matrix-free methods
14: Concepts: user-provided preconditioner;
15: Processors: 1
16: T*/
18: /*
19: Include "petscsnes.h" so that we can use SNES solvers. Note that this
20: file automatically includes:
21: petscsys.h - base PETSc routines petscvec.h - vectors
22: petscmat.h - matrices
23: petscis.h - index sets petscksp.h - Krylov subspace methods
24: petscviewer.h - viewers petscpc.h - preconditioners
25: petscksp.h - linear solvers
26: */
27: #include <petscsnes.h>
29: /*
30: User-defined routines
31: */
32: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
33: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
34: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
36: int main(int argc,char **argv)
37: {
38: SNES snes; /* SNES context */
39: KSP ksp; /* KSP context */
40: PC pc; /* PC context */
41: Vec x,r,F; /* vectors */
42: Mat J,JPrec; /* Jacobian,preconditioner matrices */
44: PetscInt it,n = 5,i;
45: PetscMPIInt size;
46: PetscInt *Shistit = 0,Khistl = 200,Shistl = 10;
47: PetscReal h,xp = 0.0,*Khist = 0,*Shist = 0;
48: PetscScalar v,pfive = .5;
49: PetscBool flg;
51: PetscInitialize(&argc,&argv,(char *)0,help);
52: MPI_Comm_size(PETSC_COMM_WORLD,&size);
53: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
54: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
55: h = 1.0/(n-1);
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create nonlinear solver context
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: SNESCreate(PETSC_COMM_WORLD,&snes);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Create vector data structures; set function evaluation routine
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: VecCreate(PETSC_COMM_SELF,&x);
68: VecSetSizes(x,PETSC_DECIDE,n);
69: VecSetFromOptions(x);
70: VecDuplicate(x,&r);
71: VecDuplicate(x,&F);
73: SNESSetFunction(snes,r,FormFunction,(void*)F);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create matrix data structures; set Jacobian evaluation routine
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&J);
80: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,1,PETSC_NULL,&JPrec);
82: /*
83: Note that in this case we create separate matrices for the Jacobian
84: and preconditioner matrix. Both of these are computed in the
85: routine FormJacobian()
86: */
87: SNESSetJacobian(snes,J,JPrec,FormJacobian,0);
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Customize nonlinear solver; set runtime options
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: /* Set preconditioner for matrix-free method */
94: flg = PETSC_FALSE;
95: PetscOptionsGetBool(PETSC_NULL,"-snes_mf",&flg,PETSC_NULL);
96: if (flg) {
97: SNESGetKSP(snes,&ksp);
98: KSPGetPC(ksp,&pc);
99: PetscOptionsHasName(PETSC_NULL,"-user_precond",&flg);
100: if (flg) { /* user-defined precond */
101: PCSetType(pc,PCSHELL);
102: PCShellSetApply(pc,MatrixFreePreconditioner);
103: } else {PCSetType(pc,PCNONE);}
104: }
106: SNESSetFromOptions(snes);
108: /*
109: Save all the linear residuals for all the Newton steps; this enables us
110: to retain complete convergence history for printing after the conclusion
111: of SNESSolve(). Alternatively, one could use the monitoring options
112: -snes_monitor -ksp_monitor
113: to see this information during the solver's execution; however, such
114: output during the run distorts performance evaluation data. So, the
115: following is a good option when monitoring code performance, for example
116: when using -log_summary.
117: */
118: PetscOptionsHasName(PETSC_NULL,"-rhistory",&flg);
119: if (flg) {
120: SNESGetKSP(snes,&ksp);
121: PetscMalloc(Khistl*sizeof(PetscReal),&Khist);
122: KSPSetResidualHistory(ksp,Khist,Khistl,PETSC_FALSE);
123: PetscMalloc(Shistl*sizeof(PetscReal),&Shist);
124: PetscMalloc(Shistl*sizeof(PetscInt),&Shistit);
125: SNESSetConvergenceHistory(snes,Shist,Shistit,Shistl,PETSC_FALSE);
126: }
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Initialize application:
130: Store right-hand-side of PDE and exact solution
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: xp = 0.0;
134: for (i=0; i<n; i++) {
135: v = 6.0*xp + pow(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
136: VecSetValues(F,1,&i,&v,INSERT_VALUES);
137: xp += h;
138: }
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Evaluate initial guess; then solve nonlinear system
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: VecSet(x,pfive);
145: SNESSolve(snes,PETSC_NULL,x);
146: SNESGetIterationNumber(snes,&it);
147: PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it);
149: PetscOptionsHasName(PETSC_NULL,"-rhistory",&flg);
150: if (flg) {
151: KSPGetResidualHistory(ksp,PETSC_NULL,&Khistl);
152: PetscRealView(Khistl,Khist,PETSC_VIEWER_STDOUT_SELF);
153: PetscFree(Khist);
154: SNESGetConvergenceHistory(snes,PETSC_NULL,PETSC_NULL,&Shistl);
155: PetscRealView(Shistl,Shist,PETSC_VIEWER_STDOUT_SELF);
156: PetscIntView(Shistl,Shistit,PETSC_VIEWER_STDOUT_SELF);
157: PetscFree(Shist);
158: PetscFree(Shistit);
159: }
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Free work space. All PETSc objects should be destroyed when they
163: are no longer needed.
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: VecDestroy(&x); VecDestroy(&r);
167: VecDestroy(&F); MatDestroy(&J);
168: MatDestroy(&JPrec); SNESDestroy(&snes);
169: PetscFinalize();
171: return 0;
172: }
173: /* ------------------------------------------------------------------- */
174: /*
175: FormInitialGuess - Forms initial approximation.
177: Input Parameters:
178: user - user-defined application context
179: X - vector
181: Output Parameter:
182: X - vector
183: */
184: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
185: {
186: PetscScalar *xx,*ff,*FF,d;
188: PetscInt i,n;
190: VecGetArray(x,&xx);
191: VecGetArray(f,&ff);
192: VecGetArray((Vec)dummy,&FF);
193: VecGetSize(x,&n);
194: d = (PetscReal)(n - 1); d = d*d;
195: ff[0] = xx[0];
196: for (i=1; i<n-1; i++) {
197: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
198: }
199: ff[n-1] = xx[n-1] - 1.0;
200: VecRestoreArray(x,&xx);
201: VecRestoreArray(f,&ff);
202: VecRestoreArray((Vec)dummy,&FF);
203: return 0;
204: }
205: /* ------------------------------------------------------------------- */
206: /*
207: FormJacobian - This routine demonstrates the use of different
208: matrices for the Jacobian and preconditioner
210: Input Parameters:
211: . snes - the SNES context
212: . x - input vector
213: . ptr - optional user-defined context, as set by SNESSetJacobian()
215: Output Parameters:
216: . A - Jacobian matrix
217: . B - different preconditioning matrix
218: . flag - flag indicating matrix structure
219: */
220: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *prejac,MatStructure *flag,void *dummy)
221: {
222: PetscScalar *xx,A[3],d;
223: PetscInt i,n,j[3];
226: VecGetArray(x,&xx);
227: VecGetSize(x,&n);
228: d = (PetscReal)(n - 1); d = d*d;
230: /* Form Jacobian. Also form a different preconditioning matrix that
231: has only the diagonal elements. */
232: i = 0; A[0] = 1.0;
233: MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
234: MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
235: for (i=1; i<n-1; i++) {
236: j[0] = i - 1; j[1] = i; j[2] = i + 1;
237: A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d;
238: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
239: MatSetValues(*prejac,1,&i,1,&i,&A[1],INSERT_VALUES);
240: }
241: i = n-1; A[0] = 1.0;
242: MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
243: MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
245: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
246: MatAssemblyBegin(*prejac,MAT_FINAL_ASSEMBLY);
247: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
248: MatAssemblyEnd(*prejac,MAT_FINAL_ASSEMBLY);
250: VecRestoreArray(x,&xx);
251: *flag = SAME_NONZERO_PATTERN;
252: return 0;
253: }
254: /* ------------------------------------------------------------------- */
255: /*
256: MatrixFreePreconditioner - This routine demonstrates the use of a
257: user-provided preconditioner. This code implements just the null
258: preconditioner, which of course is not recommended for general use.
260: Input Parameters:
261: + pc - preconditioner
262: - x - input vector
264: Output Parameter:
265: . y - preconditioned vector
266: */
267: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
268: {
270: VecCopy(x,y);
271: return 0;
272: }