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