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, (char *)0, 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_DEFAULT, PETSC_DEFAULT, 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*/