Actual source code: ex4.c
2: static char help[] = "Tests SNESLinesearch handling of Inf/Nan.\n\n";
4: /*
5: Include "petscsnes.h" so that we can use SNES solvers. Note that this
6: file automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets petscksp.h - Krylov subspace methods
10: petscviewer.h - viewers petscpc.h - preconditioners
11: petscksp.h - linear solvers
12: */
13: /*F
14: This examples solves
15: \begin{equation}
16: F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
17: \end{equation}
18: F*/
19: #include <petscsnes.h>
21: /*
22: User-defined routines
23: */
24: extern PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *);
25: extern PetscErrorCode FormFunction2(SNES, Vec, Vec, void *);
26: extern PetscErrorCode FormObjective(SNES, Vec, PetscReal *, void *);
28: /*
29: This is a very hacking way to trigger the objective function generating an infinity at a particular count to the call FormObjective().
30: Different line searches evaluate the full step at different counts. For l2 it is the third call (infatcount == 2) while for bt it is the second call.
31: */
32: PetscInt infatcount = 0;
34: int main(int argc, char **argv)
35: {
36: SNES snes; /* nonlinear solver context */
37: KSP ksp; /* linear solver context */
38: PC pc; /* preconditioner context */
39: Vec x, r; /* solution, residual vectors */
40: Mat J; /* Jacobian matrix */
41: PetscInt its;
42: PetscMPIInt size;
43: PetscScalar *xx;
44: PetscBool flg;
45: char type[256];
47: PetscFunctionBeginUser;
48: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
49: PetscCall(PetscOptionsGetString(NULL, NULL, "-snes_linesearch_type", type, sizeof(type), &flg));
50: if (flg) {
51: PetscCall(PetscStrcmp(type, SNESLINESEARCHBT, &flg));
52: if (flg) infatcount = 1;
53: PetscCall(PetscStrcmp(type, SNESLINESEARCHL2, &flg));
54: if (flg) infatcount = 2;
55: }
57: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
58: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs");
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Create nonlinear solver context
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
65: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: Create matrix and vector data structures; set corresponding routines
67: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: /*
69: Create vectors for solution and nonlinear function
70: */
71: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
72: PetscCall(VecSetSizes(x, PETSC_DECIDE, 2));
73: PetscCall(VecSetFromOptions(x));
74: PetscCall(VecDuplicate(x, &r));
76: /*
77: Create Jacobian matrix data structure
78: */
79: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
80: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
81: PetscCall(MatSetFromOptions(J));
82: PetscCall(MatSetUp(J));
84: PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL));
85: PetscCall(SNESSetObjective(snes, FormObjective, NULL));
86: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL));
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Customize nonlinear solver; set runtime options
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: /*
92: Set linear solver defaults for this problem. By extracting the
93: KSP and PC contexts from the SNES context, we can then
94: directly call any KSP and PC routines to set various options.
95: */
96: PetscCall(SNESGetKSP(snes, &ksp));
97: PetscCall(KSPGetPC(ksp, &pc));
98: PetscCall(PCSetType(pc, PCNONE));
99: PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20));
101: /*
102: Set SNES/KSP/KSP/PC runtime options, e.g.,
103: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
104: These options will override those specified above as long as
105: SNESSetFromOptions() is called _after_ any other customization
106: routines.
107: */
108: PetscCall(SNESSetFromOptions(snes));
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Evaluate initial guess; then solve nonlinear system
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: PetscCall(VecGetArray(x, &xx));
114: xx[0] = 2.0;
115: xx[1] = 3.0;
116: PetscCall(VecRestoreArray(x, &xx));
118: /*
119: Note: The user should initialize the vector, x, with the initial guess
120: for the nonlinear solver prior to calling SNESSolve(). In particular,
121: to employ an initial guess of zero, the user should explicitly set
122: this vector to zero by calling VecSet().
123: */
125: PetscCall(SNESSolve(snes, NULL, x));
126: PetscCall(SNESGetIterationNumber(snes, &its));
127: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Free work space. All PETSc objects should be destroyed when they
131: are no longer needed.
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: PetscCall(VecDestroy(&x));
135: PetscCall(VecDestroy(&r));
136: PetscCall(MatDestroy(&J));
137: PetscCall(SNESDestroy(&snes));
138: PetscCall(PetscFinalize());
139: return 0;
140: }
142: PetscErrorCode FormObjective(SNES snes, Vec x, PetscReal *f, void *dummy)
143: {
144: Vec F;
145: static PetscInt cnt = 0;
146: const PetscReal one = 1.0, zero = 0.0;
147: PetscReal inf;
149: PetscFunctionBeginUser;
150: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
151: inf = one / zero;
152: PetscCall(PetscFPTrapPop());
153: if (cnt++ == infatcount) *f = inf;
154: else {
155: PetscCall(VecDuplicate(x, &F));
156: PetscCall(FormFunction2(snes, x, F, dummy));
157: PetscCall(VecNorm(F, NORM_2, f));
158: PetscCall(VecDestroy(&F));
159: *f = (*f) * (*f);
160: }
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
165: {
166: const PetscScalar *xx;
167: PetscScalar *ff;
169: PetscFunctionBeginUser;
170: /*
171: Get pointers to vector data.
172: - For default PETSc vectors, VecGetArray() returns a pointer to
173: the data array. Otherwise, the routine is implementation dependent.
174: - You MUST call VecRestoreArray() when you no longer need access to
175: the array.
176: */
177: PetscCall(VecGetArrayRead(x, &xx));
178: PetscCall(VecGetArray(f, &ff));
180: /*
181: Compute function
182: */
183: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
184: ff[1] = xx[1];
186: /*
187: Restore vectors
188: */
189: PetscCall(VecRestoreArrayRead(x, &xx));
190: PetscCall(VecRestoreArray(f, &ff));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
195: {
196: const PetscScalar *xx;
197: PetscScalar A[4];
198: PetscInt idx[2] = {0, 1};
200: PetscFunctionBeginUser;
201: /*
202: Get pointer to vector data
203: */
204: PetscCall(VecGetArrayRead(x, &xx));
206: /*
207: Compute Jacobian entries and insert into matrix.
208: - Since this is such a small problem, we set all entries for
209: the matrix at once.
210: */
211: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
212: A[1] = 0.0;
213: A[2] = 0.0;
214: A[3] = 1.0;
215: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
217: /*
218: Restore vector
219: */
220: PetscCall(VecRestoreArrayRead(x, &xx));
222: /*
223: Assemble matrix
224: */
225: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
226: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
227: if (jac != B) {
228: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
229: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
230: }
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*TEST
236: test:
237: args: -snes_converged_reason -snes_linesearch_monitor -snes_linesearch_type l2
238: filter: grep Inf
240: TEST*/