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: /*T
  5:    Concepts: SNES^basic example
  6: T*/

  8: /*
  9:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 10:    file automatically includes:
 11:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 12:      petscmat.h - matrices
 13:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 14:      petscviewer.h - viewers               petscpc.h  - preconditioners
 15:      petscksp.h   - linear solvers
 16: */
 17: #include <petscsnes.h>

 19: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
 20: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);

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

 31:   PetscInitialize(&argc,&argv,(char*)0,help);

 33:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 34:      Create nonlinear solver context
 35:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 36:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

 49:   /*
 50:      Create Jacobian matrix data structure
 51:   */
 52:   MatCreate(PETSC_COMM_WORLD,&J);
 53:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 54:   MatSetFromOptions(J);
 55:   MatSetUp(J);

 57:   /*
 58:      Set function evaluation routine and vector.
 59:   */
 60:   SNESSetFunction(snes,r,FormFunction1,NULL);

 62:   /*
 63:      Set Jacobian matrix data structure and Jacobian evaluation routine
 64:   */
 65:   SNESSetJacobian(snes,J,J,FormJacobian1,NULL);

 67:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 68:      Customize nonlinear solver; set runtime options
 69:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 70:   SNESSetFromOptions(snes);

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

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

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

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

103:   VecDestroy(&x)); PetscCall(VecDestroy(&r);
104:   MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
105:   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:   /*
126:     Get pointers to vector data.
127:     - For default PETSc vectors, VecGetArray() returns a pointer to
128:     the data array.  Otherwise, the routine is implementation dependent.
129:     - You MUST call VecRestoreArray() when you no longer need access to
130:     the array.
131:   */
132:   VecGetArrayRead(x,&xx);
133:   VecGetArray(f,&ff);

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

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

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

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

164:   /*
165:      Get pointer to vector data
166:   */
167:   VecGetArrayRead(x,&xx);

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

180:   /*
181:      Restore vector
182:   */
183:   VecRestoreArrayRead(x,&xx);

185:   /*
186:      Assemble matrix
187:   */
188:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
189:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
190:   if (jac != B) {
191:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
192:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
193:   }
194:   return 0;
195: }

197: /*TEST

199:    test:
200:       suffix: 1
201:       args: -snes_monitor_short -snes_max_it 1000
202:       requires: !single

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

209: TEST*/