1 static const char help[] = "Attempts to solve for root of a function with multiple local minima.\n\
2 With the proper initial guess, a backtracking line-search fails because Newton's method gets\n\
3 stuck in a local minimum. However, a critical point line-search or Newton's method without a\n\
4 line search succeeds.\n";
5
6 /* Solve 1D problem f(x) = 8 * exp(-4 * (x - 2)^2) * (x - 2) + 2 * x = 0
7
8 This problem is based on the example given here: https://scicomp.stackexchange.com/a/2446/24756
9 Originally an optimization problem to find the minimum of the function
10
11 g(x) = x^2 - exp(-4 * (x - 2)^2)
12
13 it has been reformulated to solve dg(x)/dx = f(x) = 0. The reformulated problem has several local
14 minima that can cause problems for some global Newton root-finding methods. In this particular
15 example, an initial guess of x0 = 2.5 generates an initial search direction (-df/dx is positive)
16 away from the root and towards a local minimum in which a back-tracking line search gets trapped.
17 However, omitting a line-search or using a critical point or bisection line search, the solve is
18 successful.
19
20 The test outputs the final result for x and f(x).
21
22 Example usage:
23
24 Get help:
25 ./ex99 -help
26
27 Monitor run (with default back-tracking line search; solve fails):
28 ./ex99 -snes_converged_reason -snes_monitor -snes_linesearch_monitor -ksp_converged_reason -ksp_monitor
29
30 Run without a line search; solve succeeds:
31 ./ex99 -snes_linesearch_type basic
32
33 Run with a critical point line search; solve succeeds:
34 ./ex99 -snes_linesearch_type cp
35 */
36
37 #include <math.h>
38 #include <petscsnes.h>
39
40 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
41 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
42
main(int argc,char ** argv)43 int main(int argc, char **argv)
44 {
45 SNES snes; /* nonlinear solver context */
46 KSP ksp; /* linear solver context */
47 PC pc; /* preconditioner context */
48 Vec x, r; /* solution, residual vectors */
49 Mat J; /* Jacobian matrix */
50 PetscMPIInt size;
51
52 PetscFunctionBeginUser;
53 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
54 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
55 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs");
56
57 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58 Create nonlinear solver context
59 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
61
62 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 Create matrix and vector data structures; set corresponding routines
64 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 /*
66 Create vectors for solution and nonlinear function
67 */
68 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
69 PetscCall(VecSetSizes(x, PETSC_DECIDE, 1));
70 PetscCall(VecSetFromOptions(x));
71 PetscCall(VecDuplicate(x, &r));
72
73 /*
74 Create Jacobian matrix data structure
75 */
76 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
77 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
78 PetscCall(MatSetFromOptions(J));
79 PetscCall(MatSetUp(J));
80
81 PetscCall(SNESSetFunction(snes, r, FormFunction, NULL));
82 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));
83
84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85 Customize nonlinear solver; set runtime options
86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87 /*
88 Set linear solver defaults for this problem. By extracting the
89 KSP and PC contexts from the SNES context, we can then
90 directly call any KSP and PC routines to set various options.
91 */
92 PetscCall(SNESGetKSP(snes, &ksp));
93 PetscCall(KSPGetPC(ksp, &pc));
94 PetscCall(PCSetType(pc, PCNONE));
95 PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_CURRENT, PETSC_CURRENT, 20));
96
97 /*
98 Set SNES/KSP/KSP/PC runtime options, e.g.,
99 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
100 These options will override those specified above as long as
101 SNESSetFromOptions() is called _after_ any other customization
102 routines.
103 */
104 PetscCall(SNESSetFromOptions(snes));
105
106 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107 Evaluate initial guess; then solve nonlinear system
108 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109 PetscCall(VecSet(x, 2.5));
110
111 PetscCall(SNESSolve(snes, NULL, x));
112
113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114 Output x and f(x)
115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
117 PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
118
119 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120 Free work space. All PETSc objects should be destroyed when they
121 are no longer needed.
122 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123
124 PetscCall(VecDestroy(&x));
125 PetscCall(VecDestroy(&r));
126 PetscCall(MatDestroy(&J));
127 PetscCall(SNESDestroy(&snes));
128 PetscCall(PetscFinalize());
129 return 0;
130 }
131
FormFunction(SNES snes,Vec x,Vec f,PetscCtx ctx)132 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
133 {
134 const PetscScalar *xx;
135 PetscScalar *ff;
136
137 PetscFunctionBeginUser;
138 /*
139 Get pointers to vector data.
140 - For default PETSc vectors, VecGetArray() returns a pointer to
141 the data array. Otherwise, the routine is implementation dependent.
142 - You MUST call VecRestoreArray() when you no longer need access to
143 the array.
144 */
145 PetscCall(VecGetArrayRead(x, &xx));
146 PetscCall(VecGetArray(f, &ff));
147
148 /* Compute function */
149 ff[0] = 8. * PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * (xx[0] - 2.) + 2. * xx[0];
150
151 /* Restore vectors */
152 PetscCall(VecRestoreArrayRead(x, &xx));
153 PetscCall(VecRestoreArray(f, &ff));
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void * dummy)157 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
158 {
159 const PetscScalar *xx;
160 PetscScalar A[1];
161 PetscInt idx[1] = {0};
162
163 PetscFunctionBeginUser;
164 /*
165 Get pointer to vector data
166 */
167 PetscCall(VecGetArrayRead(x, &xx));
168
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] = 8. * ((xx[0] - 2.) * (PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * -8. * (xx[0] - 2.)) + PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.))) + 2.;
175
176 PetscCall(MatSetValues(B, 1, idx, 1, idx, A, INSERT_VALUES));
177
178 /*
179 Restore vector
180 */
181 PetscCall(VecRestoreArrayRead(x, &xx));
182
183 /*
184 Assemble matrix
185 */
186 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
187 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
188 if (jac != B) {
189 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
190 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
191 }
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
195 /*TEST
196
197 test:
198 suffix: 1
199 args: -snes_linesearch_type cp
200 test:
201 suffix: 2
202 args: -snes_linesearch_type basic
203 test:
204 suffix: 3
205 test:
206 suffix: 4
207 args: -snes_type newtontrdc
208 test:
209 suffix: 5
210 args: -snes_type newtontrdc -snes_trdc_use_cauchy false
211 test:
212 suffix: 6
213 args: -snes_linesearch_type bisection
214
215 TEST*/
216