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 or bisection line search, the solve is
18: successful.
20: The test outputs the final result for x and f(x).
22: Example usage:
24: Get help:
25: ./ex99 -help
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
30: Run without a line search; solve succeeds:
31: ./ex99 -snes_linesearch_type basic
33: Run with a critical point line search; solve succeeds:
34: ./ex99 -snes_linesearch_type cp
35: */
37: #include <math.h>
38: #include <petscsnes.h>
40: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
41: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
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;
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");
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create nonlinear solver context
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
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));
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));
81: PetscCall(SNESSetFunction(snes, r, FormFunction, NULL));
82: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));
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));
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));
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Evaluate initial guess; then solve nonlinear system
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: PetscCall(VecSet(x, 2.5));
111: PetscCall(SNESSolve(snes, NULL, x));
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Output x and f(x)
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
117: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Free work space. All PETSc objects should be destroyed when they
121: are no longer needed.
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: PetscCall(VecDestroy(&x));
125: PetscCall(VecDestroy(&r));
126: PetscCall(MatDestroy(&J));
127: PetscCall(SNESDestroy(&snes));
128: PetscCall(PetscFinalize());
129: return 0;
130: }
132: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
133: {
134: const PetscScalar *xx;
135: PetscScalar *ff;
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));
148: /* Compute function */
149: ff[0] = 8. * PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * (xx[0] - 2.) + 2. * xx[0];
151: /* Restore vectors */
152: PetscCall(VecRestoreArrayRead(x, &xx));
153: PetscCall(VecRestoreArray(f, &ff));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
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};
163: PetscFunctionBeginUser;
164: /*
165: Get pointer to vector data
166: */
167: PetscCall(VecGetArrayRead(x, &xx));
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.;
176: PetscCall(MatSetValues(B, 1, idx, 1, idx, A, INSERT_VALUES));
178: /*
179: Restore vector
180: */
181: PetscCall(VecRestoreArrayRead(x, &xx));
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: }
195: /*TEST
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
215: TEST*/