Actual source code: ex42.c
1: static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n";
3: /*
4: Include "petscsnes.h" so that we can use SNES solvers. Note that this
5: file automatically includes:
6: petscsys.h - base PETSc routines petscvec.h - vectors
7: petscmat.h - matrices
8: petscis.h - index sets petscksp.h - Krylov subspace methods
9: petscviewer.h - viewers petscpc.h - preconditioners
10: petscksp.h - linear solvers
11: */
12: #include <petscsnes.h>
14: extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
15: extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
17: int main(int argc, char **argv)
18: {
19: SNES snes; /* nonlinear solver context */
20: Vec x, r; /* solution, residual vectors */
21: Mat J; /* Jacobian matrix */
22: PetscInt its;
23: PetscScalar *xx;
24: SNESConvergedReason reason;
26: PetscFunctionBeginUser;
27: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
29: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: Create nonlinear solver context
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
32: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
34: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: Create matrix and vector data structures; set corresponding routines
36: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37: /*
38: Create vectors for solution and nonlinear function
39: */
40: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
41: PetscCall(VecSetSizes(x, PETSC_DECIDE, 2));
42: PetscCall(VecSetFromOptions(x));
43: PetscCall(VecDuplicate(x, &r));
45: /*
46: Create Jacobian matrix data structure
47: */
48: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
49: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
50: PetscCall(MatSetFromOptions(J));
51: PetscCall(MatSetUp(J));
53: /*
54: Set function evaluation routine and vector.
55: */
56: PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));
58: /*
59: Set Jacobian matrix data structure and Jacobian evaluation routine
60: */
61: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Customize nonlinear solver; set runtime options
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: PetscCall(SNESSetFromOptions(snes));
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Evaluate initial guess; then solve nonlinear system
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: PetscCall(VecGetArray(x, &xx));
72: xx[0] = -1.2;
73: xx[1] = 1.0;
74: PetscCall(VecRestoreArray(x, &xx));
76: /*
77: Note: The user should initialize the vector, x, with the initial guess
78: for the nonlinear solver prior to calling SNESSolve(). In particular,
79: to employ an initial guess of zero, the user should explicitly set
80: this vector to zero by calling VecSet().
81: */
83: PetscCall(SNESSolve(snes, NULL, x));
84: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
85: PetscCall(SNESGetIterationNumber(snes, &its));
86: PetscCall(SNESGetConvergedReason(snes, &reason));
87: /*
88: Some systems computes a residual that is identically zero, thus converging
89: due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
90: report the reason if the iteration did not converge so that the tests are
91: reproducible.
92: */
93: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s number of SNES iterations = %" PetscInt_FMT "\n", reason > 0 ? "CONVERGED" : SNESConvergedReasons[reason], its));
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Free work space. All PETSc objects should be destroyed when they
97: are no longer needed.
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: PetscCall(VecDestroy(&x));
101: PetscCall(VecDestroy(&r));
102: PetscCall(MatDestroy(&J));
103: PetscCall(SNESDestroy(&snes));
104: PetscCall(PetscFinalize());
105: return 0;
106: }
107: /* ------------------------------------------------------------------- */
108: /*
109: FormFunction1 - Evaluates nonlinear function, F(x).
111: Input Parameters:
112: . snes - the SNES context
113: . x - input vector
114: . ctx - optional user-defined context
116: Output Parameter:
117: . f - function vector
118: */
119: PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx)
120: {
121: PetscScalar *ff;
122: const PetscScalar *xx;
124: PetscFunctionBeginUser;
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: PetscCall(VecGetArrayRead(x, &xx));
133: PetscCall(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: PetscCall(VecRestoreArrayRead(x, &xx));
141: PetscCall(VecRestoreArray(f, &ff));
142: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBeginUser;
165: /*
166: Get pointer to vector data
167: */
168: PetscCall(VecGetArrayRead(x, &xx));
170: /*
171: Compute Jacobian entries and insert into matrix.
172: - Since this is such a small problem, we set all entries for
173: the matrix at once.
174: */
175: A[0] = 2.0 + 1200.0 * xx[0] * xx[0] - 400.0 * xx[1];
176: A[1] = -400.0 * xx[0];
177: A[2] = -400.0 * xx[0];
178: A[3] = 200;
179: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
181: /*
182: Restore vector
183: */
184: PetscCall(VecRestoreArrayRead(x, &xx));
186: /*
187: Assemble matrix
188: */
189: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
190: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
191: if (jac != B) {
192: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
193: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
194: }
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: /*TEST
200: test:
201: suffix: 1
202: args: -snes_monitor_short -snes_max_it 1000
203: requires: !single
205: test:
206: suffix: 2
207: args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false
208: requires: !single
210: TEST*/