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>


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

 34:   PetscInitialize(&argc,&argv,(char *)0,help);
 35: 
 36:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 37:      Create nonlinear solver context
 38:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 39:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

 59:   /* 
 60:      Set function evaluation routine and vector.
 61:   */
 62:   SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);

 64:   /* 
 65:      Set Jacobian matrix data structure and Jacobian evaluation routine
 66:   */
 67:   SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);

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

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

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

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

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

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

116:    Input Parameters:
117: .  snes - the SNES context
118: .  x    - input vector
119: .  ctx  - optional user-defined context

121:    Output Parameter:
122: .  f - function vector
123:  */
124: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
125: {
127:   PetscScalar    *xx,*ff;

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

139:   /* Compute function */
140:   ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
141:   ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];

143:   /* Restore vectors */
144:   VecRestoreArray(x,&xx);
145:   VecRestoreArray(f,&ff);
146:   return 0;
147: }
148: /* ------------------------------------------------------------------- */
151: /*
152:    FormJacobian1 - Evaluates Jacobian matrix.

154:    Input Parameters:
155: .  snes - the SNES context
156: .  x - input vector
157: .  dummy - optional user-defined context (not used here)

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

170:   /*
171:      Get pointer to vector data
172:   */
173:   VecGetArray(x,&xx);

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

187:   /*
188:      Restore vector
189:   */
190:   VecRestoreArray(x,&xx);

192:   /* 
193:      Assemble matrix
194:   */
195:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
196:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
197:   if (*jac != *B){
198:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
199:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
200:   }
201:   return 0;
202: }