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;

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

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

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

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

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

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

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

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

 87:   SNESSolve(snes,PETSC_NULL,x);
 88:   VecView(x,0);
 89:   SNESGetIterationNumber(snes,&its);
 90:   PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Free work space.  All PETSc objects should be destroyed when they
 94:      are no longer needed.
 95:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 97:   VecDestroy(x); VecDestroy(r);
 98:   MatDestroy(J); SNESDestroy(snes);
 99:   PetscFinalize();
100:   return 0;
101: }
102: /* ------------------------------------------------------------------- */
105: /* 
106:    FormFunction1 - Evaluates nonlinear function, F(x).

108:    Input Parameters:
109: .  snes - the SNES context
110: .  x    - input vector
111: .  ctx  - optional user-defined context

113:    Output Parameter:
114: .  f - function vector
115:  */
116: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
117: {
119:   PetscScalar    *xx,*ff;

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

131:   /* Compute function */
132:   ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
133:   ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];

135:   /* Restore vectors */
136:   VecRestoreArray(x,&xx);
137:   VecRestoreArray(f,&ff);
138:   return 0;
139: }
140: /* ------------------------------------------------------------------- */
143: /*
144:    FormJacobian1 - Evaluates Jacobian matrix.

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

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

162:   /*
163:      Get pointer to vector data
164:   */
165:   VecGetArray(x,&xx);

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

179:   /*
180:      Restore vector
181:   */
182:   VecRestoreArray(x,&xx);

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