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