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 line search, the solve is successful.
18
19 The test outputs the final result for x and f(x).
20
21 Example usage:
22
23 Get help:
24 ./ex99 -help
25
26 Monitor run (with default back-tracking line search; solve fails):
27 ./ex99 -snes_converged_reason -snes_monitor -snes_linesearch_monitor -ksp_converged_reason -ksp_monitor
28
29 Run without a line search; solve succeeds:
30 ./ex99 -snes_linesearch_type basic
31
32 Run with a critical point line search; solve succeeds:
33 ./ex99 -snes_linesearch_type cp
34 */
35
36 #include <math.h>
37 #include <petscsnes.h>
38
39 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
40 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
41
main(int argc,char ** argv)42 int main(int argc,char **argv)
43 {
44 SNES snes; /* nonlinear solver context */
45 KSP ksp; /* linear solver context */
46 PC pc; /* preconditioner context */
47 Vec x,r; /* solution, residual vectors */
48 Mat J; /* Jacobian matrix */
49 PetscErrorCode ierr;
50 PetscMPIInt size;
51
52 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
53 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
54 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");
55
56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 Create nonlinear solver context
58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
60
61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62 Create matrix and vector data structures; set corresponding routines
63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64 /*
65 Create vectors for solution and nonlinear function
66 */
67 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
68 ierr = VecSetSizes(x,PETSC_DECIDE,1);CHKERRQ(ierr);
69 ierr = VecSetFromOptions(x);CHKERRQ(ierr);
70 ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
71
72 /*
73 Create Jacobian matrix data structure
74 */
75 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
76 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,1,1);CHKERRQ(ierr);
77 ierr = MatSetFromOptions(J);CHKERRQ(ierr);
78 ierr = MatSetUp(J);CHKERRQ(ierr);
79
80 ierr = SNESSetFunction(snes,r,FormFunction,NULL);CHKERRQ(ierr);
81 ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr);
82
83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84 Customize nonlinear solver; set runtime options
85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86 /*
87 Set linear solver defaults for this problem. By extracting the
88 KSP and PC contexts from the SNES context, we can then
89 directly call any KSP and PC routines to set various options.
90 */
91 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
92 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
93 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
94 ierr = KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);CHKERRQ(ierr);
95
96 /*
97 Set SNES/KSP/KSP/PC runtime options, e.g.,
98 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
99 These options will override those specified above as long as
100 SNESSetFromOptions() is called _after_ any other customization
101 routines.
102 */
103 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
104
105 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106 Evaluate initial guess; then solve nonlinear system
107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108 ierr = VecSet(x,2.5);CHKERRQ(ierr);
109
110 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
111
112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113 Output x and f(x)
114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
116 ierr = VecView(r,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
117
118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119 Free work space. All PETSc objects should be destroyed when they
120 are no longer needed.
121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122
123 ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr);
124 ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr);
125 ierr = PetscFinalize();
126 return ierr;
127 }
128
FormFunction(SNES snes,Vec x,Vec f,void * ctx)129 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
130 {
131 PetscErrorCode ierr;
132 const PetscScalar *xx;
133 PetscScalar *ff;
134
135 /*
136 Get pointers to vector data.
137 - For default PETSc vectors, VecGetArray() returns a pointer to
138 the data array. Otherwise, the routine is implementation dependent.
139 - You MUST call VecRestoreArray() when you no longer need access to
140 the array.
141 */
142 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
143 ierr = VecGetArray(f,&ff);CHKERRQ(ierr);
144
145 /* Compute function */
146 ff[0] = 8. * PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * (xx[0] - 2.) + 2. * xx[0];
147
148 /* Restore vectors */
149 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
150 ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr);
151 return 0;
152 }
153
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void * dummy)154 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
155 {
156 const PetscScalar *xx;
157 PetscScalar A[1];
158 PetscErrorCode ierr;
159 PetscInt idx[1] = {0};
160
161 /*
162 Get pointer to vector data
163 */
164 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
165
166 /*
167 Compute Jacobian entries and insert into matrix.
168 - Since this is such a small problem, we set all entries for
169 the matrix at once.
170 */
171 A[0] = 8. * ((xx[0] - 2.) * (PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * -8. * (xx[0] - 2.))
172 + PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)))
173 + 2.;
174
175 ierr = MatSetValues(B,1,idx,1,idx,A,INSERT_VALUES);CHKERRQ(ierr);
176
177 /*
178 Restore vector
179 */
180 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
181
182 /*
183 Assemble matrix
184 */
185 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
186 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
187 if (jac != B) {
188 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
189 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
190 }
191 return 0;
192 }
193
194 /*TEST
195
196 test:
197 suffix: 1
198 args: -snes_linesearch_type cp
199 test:
200 suffix: 2
201 args: -snes_linesearch_type basic
202 test:
203 suffix: 3
204
205 TEST*/
206