Actual source code: ex42.c


  2: static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\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: #include <petscsnes.h>

 15: extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
 16: extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);

 18: int main(int argc, char **argv)
 19: {
 20:   SNES                snes; /* nonlinear solver context */
 21:   Vec                 x, r; /* solution, residual vectors */
 22:   Mat                 J;    /* Jacobian matrix */
 23:   PetscInt            its;
 24:   PetscScalar        *xx;
 25:   SNESConvergedReason reason;

 27:   PetscFunctionBeginUser;
 28:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 30:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 31:      Create nonlinear solver context
 32:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 33:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

 35:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 36:      Create matrix and vector data structures; set corresponding routines
 37:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 38:   /*
 39:      Create vectors for solution and nonlinear function
 40:   */
 41:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 42:   PetscCall(VecSetSizes(x, PETSC_DECIDE, 2));
 43:   PetscCall(VecSetFromOptions(x));
 44:   PetscCall(VecDuplicate(x, &r));

 46:   /*
 47:      Create Jacobian matrix data structure
 48:   */
 49:   PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
 50:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
 51:   PetscCall(MatSetFromOptions(J));
 52:   PetscCall(MatSetUp(J));

 54:   /*
 55:      Set function evaluation routine and vector.
 56:   */
 57:   PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));

 59:   /*
 60:      Set Jacobian matrix data structure and Jacobian evaluation routine
 61:   */
 62:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));

 64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:      Customize nonlinear solver; set runtime options
 66:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 67:   PetscCall(SNESSetFromOptions(snes));

 69:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 70:      Evaluate initial guess; then solve nonlinear system
 71:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 72:   PetscCall(VecGetArray(x, &xx));
 73:   xx[0] = -1.2;
 74:   xx[1] = 1.0;
 75:   PetscCall(VecRestoreArray(x, &xx));

 77:   /*
 78:      Note: The user should initialize the vector, x, with the initial guess
 79:      for the nonlinear solver prior to calling SNESSolve().  In particular,
 80:      to employ an initial guess of zero, the user should explicitly set
 81:      this vector to zero by calling VecSet().
 82:   */

 84:   PetscCall(SNESSolve(snes, NULL, x));
 85:   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
 86:   PetscCall(SNESGetIterationNumber(snes, &its));
 87:   PetscCall(SNESGetConvergedReason(snes, &reason));
 88:   /*
 89:      Some systems computes a residual that is identically zero, thus converging
 90:      due to FNORM_ABS, others converge due to FNORM_RELATIVE.  Here, we only
 91:      report the reason if the iteration did not converge so that the tests are
 92:      reproducible.
 93:   */
 94:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s number of SNES iterations = %" PetscInt_FMT "\n", reason > 0 ? "CONVERGED" : SNESConvergedReasons[reason], its));

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:      Free work space.  All PETSc objects should be destroyed when they
 98:      are no longer needed.
 99:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

101:   PetscCall(VecDestroy(&x));
102:   PetscCall(VecDestroy(&r));
103:   PetscCall(MatDestroy(&J));
104:   PetscCall(SNESDestroy(&snes));
105:   PetscCall(PetscFinalize());
106:   return 0;
107: }
108: /* ------------------------------------------------------------------- */
109: /*
110:    FormFunction1 - Evaluates nonlinear function, F(x).

112:    Input Parameters:
113: .  snes - the SNES context
114: .  x    - input vector
115: .  ctx  - optional user-defined context

117:    Output Parameter:
118: .  f - function vector
119:  */
120: PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx)
121: {
122:   PetscScalar       *ff;
123:   const PetscScalar *xx;

125:   PetscFunctionBeginUser;
126:   /*
127:     Get pointers to vector data.
128:     - For default PETSc vectors, VecGetArray() returns a pointer to
129:     the data array.  Otherwise, the routine is implementation dependent.
130:     - You MUST call VecRestoreArray() when you no longer need access to
131:     the array.
132:   */
133:   PetscCall(VecGetArrayRead(x, &xx));
134:   PetscCall(VecGetArray(f, &ff));

136:   /* Compute function */
137:   ff[0] = -2.0 + 2.0 * xx[0] + 400.0 * xx[0] * xx[0] * xx[0] - 400.0 * xx[0] * xx[1];
138:   ff[1] = -200.0 * xx[0] * xx[0] + 200.0 * xx[1];

140:   /* Restore vectors */
141:   PetscCall(VecRestoreArrayRead(x, &xx));
142:   PetscCall(VecRestoreArray(f, &ff));
143:   PetscFunctionReturn(PETSC_SUCCESS);
144: }
145: /* ------------------------------------------------------------------- */
146: /*
147:    FormJacobian1 - Evaluates Jacobian matrix.

149:    Input Parameters:
150: .  snes - the SNES context
151: .  x - input vector
152: .  dummy - optional user-defined context (not used here)

154:    Output Parameters:
155: .  jac - Jacobian matrix
156: .  B - optionally different preconditioning matrix
157: .  flag - flag indicating matrix structure
158: */
159: PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
160: {
161:   const PetscScalar *xx;
162:   PetscScalar        A[4];
163:   PetscInt           idx[2] = {0, 1};

165:   PetscFunctionBeginUser;
166:   /*
167:      Get pointer to vector data
168:   */
169:   PetscCall(VecGetArrayRead(x, &xx));

171:   /*
172:      Compute Jacobian entries and insert into matrix.
173:       - Since this is such a small problem, we set all entries for
174:         the matrix at once.
175:   */
176:   A[0] = 2.0 + 1200.0 * xx[0] * xx[0] - 400.0 * xx[1];
177:   A[1] = -400.0 * xx[0];
178:   A[2] = -400.0 * xx[0];
179:   A[3] = 200;
180:   PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));

182:   /*
183:      Restore vector
184:   */
185:   PetscCall(VecRestoreArrayRead(x, &xx));

187:   /*
188:      Assemble matrix
189:   */
190:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
191:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
192:   if (jac != B) {
193:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
194:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
195:   }
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: /*TEST

201:    test:
202:       suffix: 1
203:       args: -snes_monitor_short -snes_max_it 1000
204:       requires: !single

206:    test:
207:       suffix: 2
208:       args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false
209:       requires: !single

211: TEST*/