Actual source code: ex42.c
2: static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n";
4: /*T
5: Concepts: SNES^basic example
6: T*/
8: /*
9: Include "petscsnes.h" so that we can use SNES solvers. Note that this
10: file automatically includes:
11: petscsys.h - base PETSc routines petscvec.h - vectors
12: petscmat.h - matrices
13: petscis.h - index sets petscksp.h - Krylov subspace methods
14: petscviewer.h - viewers petscpc.h - preconditioners
15: petscksp.h - linear solvers
16: */
17: #include <petscsnes.h>
19: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
20: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
22: int main(int argc,char **argv)
23: {
24: SNES snes; /* nonlinear solver context */
25: Vec x,r; /* solution, residual vectors */
26: Mat J; /* Jacobian matrix */
27: PetscInt its;
28: PetscScalar *xx;
29: SNESConvergedReason reason;
31: PetscInitialize(&argc,&argv,(char*)0,help);
33: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
34: Create nonlinear solver context
35: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
36: SNESCreate(PETSC_COMM_WORLD,&snes);
38: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39: Create matrix and vector data structures; set corresponding routines
40: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41: /*
42: Create vectors for solution and nonlinear function
43: */
44: VecCreate(PETSC_COMM_WORLD,&x);
45: VecSetSizes(x,PETSC_DECIDE,2);
46: VecSetFromOptions(x);
47: VecDuplicate(x,&r);
49: /*
50: Create Jacobian matrix data structure
51: */
52: MatCreate(PETSC_COMM_WORLD,&J);
53: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
54: MatSetFromOptions(J);
55: MatSetUp(J);
57: /*
58: Set function evaluation routine and vector.
59: */
60: SNESSetFunction(snes,r,FormFunction1,NULL);
62: /*
63: Set Jacobian matrix data structure and Jacobian evaluation routine
64: */
65: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Customize nonlinear solver; set runtime options
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: SNESSetFromOptions(snes);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Evaluate initial guess; then solve nonlinear system
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: VecGetArray(x,&xx);
76: xx[0] = -1.2; xx[1] = 1.0;
77: VecRestoreArray(x,&xx);
79: /*
80: Note: The user should initialize the vector, x, with the initial guess
81: for the nonlinear solver prior to calling SNESSolve(). In particular,
82: to employ an initial guess of zero, the user should explicitly set
83: this vector to zero by calling VecSet().
84: */
86: SNESSolve(snes,NULL,x);
87: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
88: SNESGetIterationNumber(snes,&its);
89: SNESGetConvergedReason(snes,&reason);
90: /*
91: Some systems computes a residual that is identically zero, thus converging
92: due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
93: report the reason if the iteration did not converge so that the tests are
94: reproducible.
95: */
96: PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],its);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Free work space. All PETSc objects should be destroyed when they
100: are no longer needed.
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: VecDestroy(&x)); PetscCall(VecDestroy(&r);
104: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
105: PetscFinalize();
106: return 0;
107: }
108: /* ------------------------------------------------------------------- */
109: /*
110: FormFunction1 - Evaluates nonlinear function, F(x).
112: Input Parameters:
113: . snes - the SNES context
114: . x - input vector
115: . ctx - optional user-defined context
117: Output Parameter:
118: . f - function vector
119: */
120: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
121: {
122: PetscScalar *ff;
123: const PetscScalar *xx;
125: /*
126: Get pointers to vector data.
127: - For default PETSc vectors, VecGetArray() returns a pointer to
128: the data array. Otherwise, the routine is implementation dependent.
129: - You MUST call VecRestoreArray() when you no longer need access to
130: the array.
131: */
132: VecGetArrayRead(x,&xx);
133: VecGetArray(f,&ff);
135: /* Compute function */
136: ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
137: ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];
139: /* Restore vectors */
140: VecRestoreArrayRead(x,&xx);
141: VecRestoreArray(f,&ff);
142: return 0;
143: }
144: /* ------------------------------------------------------------------- */
145: /*
146: FormJacobian1 - Evaluates Jacobian matrix.
148: Input Parameters:
149: . snes - the SNES context
150: . x - input vector
151: . dummy - optional user-defined context (not used here)
153: Output Parameters:
154: . jac - Jacobian matrix
155: . B - optionally different preconditioning matrix
156: . flag - flag indicating matrix structure
157: */
158: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
159: {
160: const PetscScalar *xx;
161: PetscScalar A[4];
162: PetscInt idx[2] = {0,1};
164: /*
165: Get pointer to vector data
166: */
167: VecGetArrayRead(x,&xx);
169: /*
170: Compute Jacobian entries and insert into matrix.
171: - Since this is such a small problem, we set all entries for
172: the matrix at once.
173: */
174: A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1];
175: A[1] = -400.0*xx[0];
176: A[2] = -400.0*xx[0];
177: A[3] = 200;
178: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
180: /*
181: Restore vector
182: */
183: VecRestoreArrayRead(x,&xx);
185: /*
186: Assemble matrix
187: */
188: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
189: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
190: if (jac != B) {
191: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
192: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
193: }
194: return 0;
195: }
197: /*TEST
199: test:
200: suffix: 1
201: args: -snes_monitor_short -snes_max_it 1000
202: requires: !single
204: test:
205: suffix: 2
206: args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false
207: requires: !single
209: TEST*/