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;
51: PetscFunctionBeginUser;
52: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
53: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
54: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs");
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Create nonlinear solver context
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: PetscCall(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: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
68: PetscCall(VecSetSizes(x, PETSC_DECIDE, 1));
69: PetscCall(VecSetFromOptions(x));
70: PetscCall(VecDuplicate(x, &r));
72: /*
73: Create Jacobian matrix data structure
74: */
75: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
76: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
77: PetscCall(MatSetFromOptions(J));
78: PetscCall(MatSetUp(J));
80: PetscCall(SNESSetFunction(snes, r, FormFunction, NULL));
81: PetscCall(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: PetscCall(SNESGetKSP(snes, &ksp));
92: PetscCall(KSPGetPC(ksp, &pc));
93: PetscCall(PCSetType(pc, PCNONE));
94: PetscCall(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: PetscCall(SNESSetFromOptions(snes));
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Evaluate initial guess; then solve nonlinear system
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: PetscCall(VecSet(x, 2.5));
110: PetscCall(SNESSolve(snes, NULL, x));
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Output x and f(x)
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
116: PetscCall(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: PetscCall(VecDestroy(&x));
124: PetscCall(VecDestroy(&r));
125: PetscCall(MatDestroy(&J));
126: PetscCall(SNESDestroy(&snes));
127: PetscCall(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: PetscFunctionBeginUser;
137: /*
138: Get pointers to vector data.
139: - For default PETSc vectors, VecGetArray() returns a pointer to
140: the data array. Otherwise, the routine is implementation dependent.
141: - You MUST call VecRestoreArray() when you no longer need access to
142: the array.
143: */
144: PetscCall(VecGetArrayRead(x, &xx));
145: PetscCall(VecGetArray(f, &ff));
147: /* Compute function */
148: ff[0] = 8. * PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * (xx[0] - 2.) + 2. * xx[0];
150: /* Restore vectors */
151: PetscCall(VecRestoreArrayRead(x, &xx));
152: PetscCall(VecRestoreArray(f, &ff));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
157: {
158: const PetscScalar *xx;
159: PetscScalar A[1];
160: PetscInt idx[1] = {0};
162: PetscFunctionBeginUser;
163: /*
164: Get pointer to vector data
165: */
166: PetscCall(VecGetArrayRead(x, &xx));
168: /*
169: Compute Jacobian entries and insert into matrix.
170: - Since this is such a small problem, we set all entries for
171: the matrix at once.
172: */
173: 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: PetscCall(MatSetValues(B, 1, idx, 1, idx, A, INSERT_VALUES));
177: /*
178: Restore vector
179: */
180: PetscCall(VecRestoreArrayRead(x, &xx));
182: /*
183: Assemble matrix
184: */
185: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
186: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
187: if (jac != B) {
188: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
189: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
190: }
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*TEST
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: test:
205: suffix: 4
206: args: -snes_type newtontrdc
207: test:
208: suffix: 5
209: args: -snes_type newtontrdc -snes_trdc_use_cauchy false
211: TEST*/