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:    */
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 */
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:   */

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:   */

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*/
```