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;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
30: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: Create nonlinear solver context
32: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: PetscCall(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: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
42: PetscCall(VecSetSizes(x, PETSC_DECIDE, 2));
43: PetscCall(VecSetFromOptions(x));
44: PetscCall(VecDuplicate(x, &r));
46: /*
47: Create Jacobian matrix data structure
48: */
49: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
50: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
51: PetscCall(MatSetFromOptions(J));
52: PetscCall(MatSetUp(J));
54: /*
55: Set function evaluation routine and vector.
56: */
57: PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));
59: /*
60: Set Jacobian matrix data structure and Jacobian evaluation routine
61: */
62: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
64: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: Customize nonlinear solver; set runtime options
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: PetscCall(SNESSetFromOptions(snes));
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Evaluate initial guess; then solve nonlinear system
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: PetscCall(VecGetArray(x, &xx));
73: xx[0] = -1.2;
74: xx[1] = 1.0;
75: PetscCall(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: PetscCall(SNESSolve(snes, NULL, x));
85: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
86: PetscCall(SNESGetIterationNumber(snes, &its));
87: PetscCall(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: PetscCall(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: PetscCall(VecDestroy(&x));
102: PetscCall(VecDestroy(&r));
103: PetscCall(MatDestroy(&J));
104: PetscCall(SNESDestroy(&snes));
105: PetscCall(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: PetscFunctionBeginUser;
126: /*
127: Get pointers to vector data.
128: - For default PETSc vectors, VecGetArray() returns a pointer to
129: the data array. Otherwise, the routine is implementation dependent.
130: - You MUST call VecRestoreArray() when you no longer need access to
131: the array.
132: */
133: PetscCall(VecGetArrayRead(x, &xx));
134: PetscCall(VecGetArray(f, &ff));
136: /* Compute function */
137: ff[0] = -2.0 + 2.0 * xx[0] + 400.0 * xx[0] * xx[0] * xx[0] - 400.0 * xx[0] * xx[1];
138: ff[1] = -200.0 * xx[0] * xx[0] + 200.0 * xx[1];
140: /* Restore vectors */
141: PetscCall(VecRestoreArrayRead(x, &xx));
142: PetscCall(VecRestoreArray(f, &ff));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
145: /* ------------------------------------------------------------------- */
146: /*
147: FormJacobian1 - Evaluates Jacobian matrix.
149: Input Parameters:
150: . snes - the SNES context
151: . x - input vector
152: . dummy - optional user-defined context (not used here)
154: Output Parameters:
155: . jac - Jacobian matrix
156: . B - optionally different preconditioning matrix
157: . flag - flag indicating matrix structure
158: */
159: PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
160: {
161: const PetscScalar *xx;
162: PetscScalar A[4];
163: PetscInt idx[2] = {0, 1};
165: PetscFunctionBeginUser;
166: /*
167: Get pointer to vector data
168: */
169: PetscCall(VecGetArrayRead(x, &xx));
171: /*
172: Compute Jacobian entries and insert into matrix.
173: - Since this is such a small problem, we set all entries for
174: the matrix at once.
175: */
176: A[0] = 2.0 + 1200.0 * xx[0] * xx[0] - 400.0 * xx[1];
177: A[1] = -400.0 * xx[0];
178: A[2] = -400.0 * xx[0];
179: A[3] = 200;
180: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
182: /*
183: Restore vector
184: */
185: PetscCall(VecRestoreArrayRead(x, &xx));
187: /*
188: Assemble matrix
189: */
190: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
191: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
192: if (jac != B) {
193: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
194: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
195: }
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*TEST
201: test:
202: suffix: 1
203: args: -snes_monitor_short -snes_max_it 1000
204: requires: !single
206: test:
207: suffix: 2
208: args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false
209: requires: !single
211: TEST*/