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