Actual source code: ex4.c
1: static char help[] = "Tests SNESLinesearch handling of Inf/Nan.\n\n";
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>
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 *);
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;
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];
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, SNESLINESEARCHL2, &flg));
53: if (flg) infatcount = 2;
54: }
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");
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Create nonlinear solver context
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
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));
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));
83: PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL));
84: PetscCall(SNESSetObjective(snes, FormObjective, NULL));
85: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL));
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));
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));
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));
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: */
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));
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Free work space. All PETSc objects should be destroyed when they
130: are no longer needed.
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: PetscCall(VecDestroy(&x));
134: PetscCall(VecDestroy(&r));
135: PetscCall(MatDestroy(&J));
136: PetscCall(SNESDestroy(&snes));
137: PetscCall(PetscFinalize());
138: return 0;
139: }
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;
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: }
163: PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
164: {
165: const PetscScalar *xx;
166: PetscScalar *ff;
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));
179: /*
180: Compute function
181: */
182: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
183: ff[1] = xx[1];
185: /*
186: Restore vectors
187: */
188: PetscCall(VecRestoreArrayRead(x, &xx));
189: PetscCall(VecRestoreArray(f, &ff));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
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};
199: PetscFunctionBeginUser;
200: /*
201: Get pointer to vector data
202: */
203: PetscCall(VecGetArrayRead(x, &xx));
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));
216: /*
217: Restore vector
218: */
219: PetscCall(VecRestoreArrayRead(x, &xx));
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: }
233: /*TEST
235: test:
236: args: -snes_converged_reason -snes_linesearch_monitor -snes_linesearch_type l2
237: filter: grep Inf
239: TEST*/