Actual source code: ex99.c
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";
6: /* Solve 1D problem f(x) = 8 * exp(-4 * (x - 2)^2) * (x - 2) + 2 * x = 0
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
11: g(x) = x^2 - exp(-4 * (x - 2)^2)
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.
19: The test outputs the final result for x and f(x).
21: Example usage:
23: Get help:
24: ./ex99 -help
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
29: Run without a line search; solve succeeds:
30: ./ex99 -snes_linesearch_type basic
32: Run with a critical point line search; solve succeeds:
33: ./ex99 -snes_linesearch_type cp
34: */
36: #include <math.h>
37: #include <petscsnes.h>
39: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
40: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
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: PetscMPIInt size;
52: PetscInitialize(&argc, &argv, (char *)0, help);
53: MPI_Comm_size(PETSC_COMM_WORLD, &size);
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Create nonlinear solver context
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: SNESCreate(PETSC_COMM_WORLD, &snes);
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Create matrix and vector data structures; set corresponding routines
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: /*
65: Create vectors for solution and nonlinear function
66: */
67: VecCreate(PETSC_COMM_WORLD, &x);
68: VecSetSizes(x, PETSC_DECIDE, 1);
69: VecSetFromOptions(x);
70: VecDuplicate(x, &r);
72: /*
73: Create Jacobian matrix data structure
74: */
75: MatCreate(PETSC_COMM_WORLD, &J);
76: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
77: MatSetFromOptions(J);
78: MatSetUp(J);
80: SNESSetFunction(snes, r, FormFunction, NULL);
81: SNESSetJacobian(snes, J, J, FormJacobian, NULL);
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: SNESGetKSP(snes, &ksp);
92: KSPGetPC(ksp, &pc);
93: PCSetType(pc, PCNONE);
94: KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20);
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: SNESSetFromOptions(snes);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Evaluate initial guess; then solve nonlinear system
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: VecSet(x, 2.5);
110: SNESSolve(snes, NULL, x);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Output x and f(x)
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
116: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Free work space. All PETSc objects should be destroyed when they
120: are no longer needed.
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: VecDestroy(&x);
124: VecDestroy(&r);
125: MatDestroy(&J);
126: SNESDestroy(&snes);
127: PetscFinalize();
128: return 0;
129: }
131: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
132: {
133: const PetscScalar *xx;
134: PetscScalar *ff;
136: /*
137: Get pointers to vector data.
138: - For default PETSc vectors, VecGetArray() returns a pointer to
139: the data array. Otherwise, the routine is implementation dependent.
140: - You MUST call VecRestoreArray() when you no longer need access to
141: the array.
142: */
143: VecGetArrayRead(x, &xx);
144: VecGetArray(f, &ff);
146: /* Compute function */
147: ff[0] = 8. * PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * (xx[0] - 2.) + 2. * xx[0];
149: /* Restore vectors */
150: VecRestoreArrayRead(x, &xx);
151: VecRestoreArray(f, &ff);
152: return 0;
153: }
155: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
156: {
157: const PetscScalar *xx;
158: PetscScalar A[1];
159: PetscInt idx[1] = {0};
161: /*
162: Get pointer to vector data
163: */
164: VecGetArrayRead(x, &xx);
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.)) + PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.))) + 2.;
173: MatSetValues(B, 1, idx, 1, idx, A, INSERT_VALUES);
175: /*
176: Restore vector
177: */
178: VecRestoreArrayRead(x, &xx);
180: /*
181: Assemble matrix
182: */
183: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
184: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
185: if (jac != B) {
186: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
187: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
188: }
189: return 0;
190: }
192: /*TEST
194: test:
195: suffix: 1
196: args: -snes_linesearch_type cp
197: test:
198: suffix: 2
199: args: -snes_linesearch_type basic
200: test:
201: suffix: 3
202: test:
203: suffix: 4
204: args: -snes_type newtontrdc
205: test:
206: suffix: 5
207: args: -snes_type newtontrdc -snes_trdc_use_cauchy false
209: TEST*/