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: /*
5: Include "petscsnes.h" so that we can use SNES solvers. Note that this
6: file automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets petscksp.h - Krylov subspace methods
10: petscviewer.h - viewers petscpc.h - preconditioners
11: petscksp.h - linear solvers
12: */
13: #include <petscsnes.h>
15: extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
16: extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
18: int main(int argc, char **argv)
19: {
20: SNES snes; /* nonlinear solver context */
21: Vec x, r; /* solution, residual vectors */
22: Mat J; /* Jacobian matrix */
23: PetscInt its;
24: PetscScalar *xx;
25: SNESConvergedReason reason;
28: PetscInitialize(&argc, &argv, (char *)0, help);
30: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: Create nonlinear solver context
32: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: SNESCreate(PETSC_COMM_WORLD, &snes);
35: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: Create matrix and vector data structures; set corresponding routines
37: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38: /*
39: Create vectors for solution and nonlinear function
40: */
41: VecCreate(PETSC_COMM_WORLD, &x);
42: VecSetSizes(x, PETSC_DECIDE, 2);
43: VecSetFromOptions(x);
44: VecDuplicate(x, &r);
46: /*
47: Create Jacobian matrix data structure
48: */
49: MatCreate(PETSC_COMM_WORLD, &J);
50: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2);
51: MatSetFromOptions(J);
52: MatSetUp(J);
54: /*
55: Set function evaluation routine and vector.
56: */
57: SNESSetFunction(snes, r, FormFunction1, NULL);
59: /*
60: Set Jacobian matrix data structure and Jacobian evaluation routine
61: */
62: SNESSetJacobian(snes, J, J, FormJacobian1, NULL);
64: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: Customize nonlinear solver; set runtime options
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: SNESSetFromOptions(snes);
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Evaluate initial guess; then solve nonlinear system
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: VecGetArray(x, &xx);
73: xx[0] = -1.2;
74: xx[1] = 1.0;
75: VecRestoreArray(x, &xx);
77: /*
78: Note: The user should initialize the vector, x, with the initial guess
79: for the nonlinear solver prior to calling SNESSolve(). In particular,
80: to employ an initial guess of zero, the user should explicitly set
81: this vector to zero by calling VecSet().
82: */
84: SNESSolve(snes, NULL, x);
85: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
86: SNESGetIterationNumber(snes, &its);
87: SNESGetConvergedReason(snes, &reason);
88: /*
89: Some systems computes a residual that is identically zero, thus converging
90: due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
91: report the reason if the iteration did not converge so that the tests are
92: reproducible.
93: */
94: PetscPrintf(PETSC_COMM_WORLD, "%s number of SNES iterations = %" PetscInt_FMT "\n", reason > 0 ? "CONVERGED" : SNESConvergedReasons[reason], its);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Free work space. All PETSc objects should be destroyed when they
98: are no longer needed.
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: VecDestroy(&x);
102: VecDestroy(&r);
103: MatDestroy(&J);
104: 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*/