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*/