Actual source code: ex4.c
2: static char help[] = "Tests TSLINESEARCHL2 handing 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 either
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];
48: PetscInitialize(&argc, &argv, (char *)0, help);
49: PetscOptionsGetString(NULL, NULL, "-snes_linesearch_type", type, sizeof(type), &flg);
50: if (flg) {
51: PetscStrcmp(type, SNESLINESEARCHBT, &flg);
52: if (flg) infatcount = 1;
53: PetscStrcmp(type, SNESLINESEARCHL2, &flg);
54: if (flg) infatcount = 2;
55: }
57: MPI_Comm_size(PETSC_COMM_WORLD, &size);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Create nonlinear solver context
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: 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: VecCreate(PETSC_COMM_WORLD, &x);
72: VecSetSizes(x, PETSC_DECIDE, 2);
73: VecSetFromOptions(x);
74: VecDuplicate(x, &r);
76: /*
77: Create Jacobian matrix data structure
78: */
79: MatCreate(PETSC_COMM_WORLD, &J);
80: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2);
81: MatSetFromOptions(J);
82: MatSetUp(J);
84: SNESSetFunction(snes, r, FormFunction2, NULL);
85: SNESSetObjective(snes, FormObjective, NULL);
86: 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: SNESGetKSP(snes, &ksp);
97: KSPGetPC(ksp, &pc);
98: PCSetType(pc, PCNONE);
99: 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: SNESSetFromOptions(snes);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Evaluate initial guess; then solve nonlinear system
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: VecGetArray(x, &xx);
114: xx[0] = 2.0;
115: xx[1] = 3.0;
116: 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: SNESSolve(snes, NULL, x);
126: SNESGetIterationNumber(snes, &its);
127: 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: VecDestroy(&x);
135: VecDestroy(&r);
136: MatDestroy(&J);
137: SNESDestroy(&snes);
138: 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;
147: if (cnt++ == infatcount) *f = INFINITY;
148: else {
149: VecDuplicate(x, &F);
150: FormFunction2(snes, x, F, dummy);
151: VecNorm(F, NORM_2, f);
152: VecDestroy(&F);
153: *f = (*f) * (*f);
154: }
155: return 0;
156: }
158: PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
159: {
160: const PetscScalar *xx;
161: PetscScalar *ff;
163: /*
164: Get pointers to vector data.
165: - For default PETSc vectors, VecGetArray() returns a pointer to
166: the data array. Otherwise, the routine is implementation dependent.
167: - You MUST call VecRestoreArray() when you no longer need access to
168: the array.
169: */
170: VecGetArrayRead(x, &xx);
171: VecGetArray(f, &ff);
173: /*
174: Compute function
175: */
176: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
177: ff[1] = xx[1];
179: /*
180: Restore vectors
181: */
182: VecRestoreArrayRead(x, &xx);
183: VecRestoreArray(f, &ff);
184: return 0;
185: }
187: PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
188: {
189: const PetscScalar *xx;
190: PetscScalar A[4];
191: PetscInt idx[2] = {0, 1};
193: /*
194: Get pointer to vector data
195: */
196: VecGetArrayRead(x, &xx);
198: /*
199: Compute Jacobian entries and insert into matrix.
200: - Since this is such a small problem, we set all entries for
201: the matrix at once.
202: */
203: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
204: A[1] = 0.0;
205: A[2] = 0.0;
206: A[3] = 1.0;
207: MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES);
209: /*
210: Restore vector
211: */
212: VecRestoreArrayRead(x, &xx);
214: /*
215: Assemble matrix
216: */
217: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
218: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
219: if (jac != B) {
220: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
221: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
222: }
223: return 0;
224: }
226: /*TEST
228: build:
229: requires: infinity
231: test:
232: args: -snes_converged_reason -snes_linesearch_monitor -snes_linesearch_type l2
233: filter: grep Inf
235: TEST*/