Actual source code: ex1.c
1: /*$Id: ex1.c,v 1.26 2001/08/07 03:04:16 balay Exp $*/
3: static char help[] = "Newton's method to solve a two-variable system, sequentially.\n\n";
5: /*T
6: Concepts: SNES^basic uniprocessor example
7: Processors: 1
8: T*/
10: /*
11: Include "petscsnes.h" so that we can use SNES solvers. Note that this
12: file automatically includes:
13: petsc.h - base PETSc routines petscvec.h - vectors
14: petscsys.h - system routines petscmat.h - matrices
15: petscis.h - index sets petscksp.h - Krylov subspace methods
16: petscviewer.h - viewers petscpc.h - preconditioners
17: petscksp.h - linear solvers
18: */
19: #include petscsnes.h
21: /*
22: User-defined routines
23: */
24: extern int FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
25: extern int FormFunction1(SNES,Vec,Vec,void*);
26: extern int FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
27: extern int FormFunction2(SNES,Vec,Vec,void*);
31: int main(int argc,char **argv)
32: {
33: SNES snes; /* nonlinear solver context */
34: KSP ksp; /* linear solver context */
35: PC pc; /* preconditioner context */
36: Vec x,r; /* solution, residual vectors */
37: Mat J; /* Jacobian matrix */
38: int ierr,its,size;
39: PetscScalar pfive = .5,*xx;
40: PetscTruth flg;
42: PetscInitialize(&argc,&argv,(char *)0,help);
43: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Create nonlinear solver context
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: SNESCreate(PETSC_COMM_WORLD,&snes);
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Create matrix and vector data structures; set corresponding routines
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: /*
57: Create vectors for solution and nonlinear function
58: */
59: VecCreateSeq(PETSC_COMM_SELF,2,&x);
60: VecDuplicate(x,&r);
62: /*
63: Create Jacobian matrix data structure
64: */
65: MatCreate(PETSC_COMM_SELF,PETSC_DECIDE,PETSC_DECIDE,2,2,&J);
66: MatSetFromOptions(J);
68: PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
69: if (!flg) {
70: /*
71: Set function evaluation routine and vector.
72: */
73: SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);
75: /*
76: Set Jacobian matrix data structure and Jacobian evaluation routine
77: */
78: SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
79: } else {
80: SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);
81: SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);
82: }
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Customize nonlinear solver; set runtime options
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: /*
89: Set linear solver defaults for this problem. By extracting the
90: KSP, KSP, and PC contexts from the SNES context, we can then
91: directly call any KSP, KSP, and PC routines to set various options.
92: */
93: SNESGetKSP(snes,&ksp);
94: KSPGetPC(ksp,&pc);
95: PCSetType(pc,PCNONE);
96: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
98: /*
99: Set SNES/KSP/KSP/PC runtime options, e.g.,
100: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
101: These options will override those specified above as long as
102: SNESSetFromOptions() is called _after_ any other customization
103: routines.
104: */
105: SNESSetFromOptions(snes);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Evaluate initial guess; then solve nonlinear system
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: if (!flg) {
111: VecSet(&pfive,x);
112: } else {
113: VecGetArray(x,&xx);
114: xx[0] = 2.0; xx[1] = 3.0;
115: VecRestoreArray(x,&xx);
116: }
117: /*
118: Note: The user should initialize the vector, x, with the initial guess
119: for the nonlinear solver prior to calling SNESSolve(). In particular,
120: to employ an initial guess of zero, the user should explicitly set
121: this vector to zero by calling VecSet().
122: */
124: SNESSolve(snes,x);
125: SNESGetIterationNumber(snes,&its);
126: if (flg) {
127: Vec f;
128: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
129: SNESGetFunction(snes,&f,0,0);
130: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
131: }
133: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %d\n\n",its);
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Free work space. All PETSc objects should be destroyed when they
137: are no longer needed.
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: VecDestroy(x); VecDestroy(r);
141: MatDestroy(J); SNESDestroy(snes);
143: PetscFinalize();
144: return 0;
145: }
146: /* ------------------------------------------------------------------- */
149: /*
150: FormFunction1 - Evaluates nonlinear function, F(x).
152: Input Parameters:
153: . snes - the SNES context
154: . x - input vector
155: . dummy - optional user-defined context (not used here)
157: Output Parameter:
158: . f - function vector
159: */
160: int FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
161: {
162: int ierr;
163: PetscScalar *xx,*ff;
165: /*
166: Get pointers to vector data.
167: - For default PETSc vectors, VecGetArray() returns a pointer to
168: the data array. Otherwise, the routine is implementation dependent.
169: - You MUST call VecRestoreArray() when you no longer need access to
170: the array.
171: */
172: VecGetArray(x,&xx);
173: VecGetArray(f,&ff);
175: /*
176: Compute function
177: */
178: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
179: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
181: /*
182: Restore vectors
183: */
184: VecRestoreArray(x,&xx);
185: VecRestoreArray(f,&ff);
187: return 0;
188: }
189: /* ------------------------------------------------------------------- */
192: /*
193: FormJacobian1 - Evaluates Jacobian matrix.
195: Input Parameters:
196: . snes - the SNES context
197: . x - input vector
198: . dummy - optional user-defined context (not used here)
200: Output Parameters:
201: . jac - Jacobian matrix
202: . B - optionally different preconditioning matrix
203: . flag - flag indicating matrix structure
204: */
205: int FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
206: {
207: PetscScalar *xx,A[4];
208: int ierr,idx[2] = {0,1};
210: /*
211: Get pointer to vector data
212: */
213: VecGetArray(x,&xx);
215: /*
216: Compute Jacobian entries and insert into matrix.
217: - Since this is such a small problem, we set all entries for
218: the matrix at once.
219: */
220: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
221: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
222: MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
223: *flag = SAME_NONZERO_PATTERN;
225: /*
226: Restore vector
227: */
228: VecRestoreArray(x,&xx);
230: /*
231: Assemble matrix
232: */
233: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
234: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
236: return 0;
237: }
240: /* ------------------------------------------------------------------- */
243: int FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
244: {
245: int ierr;
246: PetscScalar *xx,*ff;
248: /*
249: Get pointers to vector data.
250: - For default PETSc vectors, VecGetArray() returns a pointer to
251: the data array. Otherwise, the routine is implementation dependent.
252: - You MUST call VecRestoreArray() when you no longer need access to
253: the array.
254: */
255: VecGetArray(x,&xx);
256: VecGetArray(f,&ff);
258: /*
259: Compute function
260: */
261: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
262: ff[1] = xx[1];
264: /*
265: Restore vectors
266: */
267: VecRestoreArray(x,&xx);
268: VecRestoreArray(f,&ff);
270: return 0;
271: }
272: /* ------------------------------------------------------------------- */
275: int FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
276: {
277: PetscScalar *xx,A[4];
278: int ierr,idx[2] = {0,1};
280: /*
281: Get pointer to vector data
282: */
283: VecGetArray(x,&xx);
285: /*
286: Compute Jacobian entries and insert into matrix.
287: - Since this is such a small problem, we set all entries for
288: the matrix at once.
289: */
290: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
291: A[2] = 0.0; A[3] = 1.0;
292: MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
293: *flag = SAME_NONZERO_PATTERN;
295: /*
296: Restore vector
297: */
298: VecRestoreArray(x,&xx);
300: /*
301: Assemble matrix
302: */
303: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
304: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
306: return 0;
307: }