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