2: static char help[] = "Newton's method for a two-variable system, sequential.\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: */

 27:  #include <petscsnes.h>

 29: /*
 30:    User-defined routines
 31: */
 32: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
 33: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 34: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
 35: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);

 37: int main(int argc,char **argv)
 38: {
 39:   SNES           snes;         /* nonlinear solver context */
 40:   KSP            ksp;          /* linear solver context */
 41:   PC             pc;           /* preconditioner context */
 42:   Vec            x,r;          /* solution, residual vectors */
 43:   Mat            J;            /* Jacobian matrix */
 45:   PetscInt       its;
 46:   PetscMPIInt    size;
 47:   PetscScalar    pfive = .5,*xx;
 48:   PetscBool      flg;

 50:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 51:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 52:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:      Create nonlinear solver context
 56:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 57:   SNESCreate(PETSC_COMM_WORLD,&snes);

 59:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 60:      Create matrix and vector data structures; set corresponding routines
 61:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 62:   /*
 63:      Create vectors for solution and nonlinear function
 64:   */
 65:   VecCreate(PETSC_COMM_WORLD,&x);
 66:   VecSetSizes(x,PETSC_DECIDE,2);
 67:   VecSetFromOptions(x);
 68:   VecDuplicate(x,&r);

 70:   /*
 71:      Create Jacobian matrix data structure
 72:   */
 73:   MatCreate(PETSC_COMM_WORLD,&J);
 74:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 75:   MatSetFromOptions(J);
 76:   MatSetUp(J);

 78:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
 79:   if (!flg) {
 80:     /*
 81:      Set function evaluation routine and vector.
 82:     */
 83:     SNESSetFunction(snes,r,FormFunction1,NULL);

 85:     /*
 86:      Set Jacobian matrix data structure and Jacobian evaluation routine
 87:     */
 88:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
 89:   } else {
 90:     SNESSetFunction(snes,r,FormFunction2,NULL);
 91:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
 92:   }

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:      Customize nonlinear solver; set runtime options
 96:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 97:   /*
 98:      Set linear solver defaults for this problem. By extracting the
 99:      KSP and PC contexts from the SNES context, we can then
100:      directly call any KSP and PC routines to set various options.
101:   */
102:   SNESGetKSP(snes,&ksp);
103:   KSPGetPC(ksp,&pc);
104:   PCSetType(pc,PCNONE);
105:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

107:   /*
108:      Set SNES/KSP/KSP/PC runtime options, e.g.,
109:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
110:      These options will override those specified above as long as
111:      SNESSetFromOptions() is called _after_ any other customization
112:      routines.
113:   */
114:   SNESSetFromOptions(snes);

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Evaluate initial guess; then solve nonlinear system
118:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119:   if (!flg) {
120:     VecSet(x,pfive);
121:   } else {
122:     VecGetArray(x,&xx);
123:     xx[0] = 2.0; xx[1] = 3.0;
124:     VecRestoreArray(x,&xx);
125:   }
126:   /*
127:      Note: The user should initialize the vector, x, with the initial guess
128:      for the nonlinear solver prior to calling SNESSolve().  In particular,
129:      to employ an initial guess of zero, the user should explicitly set
130:      this vector to zero by calling VecSet().
131:   */

133:   SNESSolve(snes,NULL,x);
134:   SNESGetIterationNumber(snes,&its);
135:   if (flg) {
136:     Vec f;
137:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
138:     SNESGetFunction(snes,&f,0,0);
139:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
140:   }

142:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:      Free work space.  All PETSc objects should be destroyed when they
146:      are no longer needed.
147:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

149:   VecDestroy(&x); VecDestroy(&r);
150:   MatDestroy(&J); SNESDestroy(&snes);
151:   PetscFinalize();
152:   return ierr;
153: }
154: /* ------------------------------------------------------------------- */
155: /*
156:    FormFunction1 - Evaluates nonlinear function, F(x).

158:    Input Parameters:
159: .  snes - the SNES context
160: .  x    - input vector
161: .  ctx  - optional user-defined context

163:    Output Parameter:
164: .  f - function vector
165:  */
166: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
167: {
168:   PetscErrorCode    ierr;
169:   const PetscScalar *xx;
170:   PetscScalar       *ff;

172:   /*
173:    Get pointers to vector data.
174:       - For default PETSc vectors, VecGetArray() returns a pointer to
175:         the data array.  Otherwise, the routine is implementation dependent.
176:       - You MUST call VecRestoreArray() when you no longer need access to
177:         the array.
178:    */
179:   VecGetArrayRead(x,&xx);
180:   VecGetArray(f,&ff);

182:   /* Compute function */
183:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
184:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

186:   /* Restore vectors */
187:   VecRestoreArrayRead(x,&xx);
188:   VecRestoreArray(f,&ff);
189:   return 0;
190: }
191: /* ------------------------------------------------------------------- */
192: /*
193:    FormJacobian1 - Evaluates Jacobian matrix.

195:    Input Parameters:
196: .  snes - the SNES context
197: .  x - input vector
198: .  dummy - optional user-defined context (not used here)

200:    Output Parameters:
201: .  jac - Jacobian matrix
202: .  B - optionally different preconditioning matrix
203: .  flag - flag indicating matrix structure
204: */
205: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
206: {
207:   const PetscScalar *xx;
208:   PetscScalar       A[4];
209:   PetscErrorCode    ierr;
210:   PetscInt          idx[2] = {0,1};

212:   /*
213:      Get pointer to vector data
214:   */
215:   VecGetArrayRead(x,&xx);

217:   /*
218:      Compute Jacobian entries and insert into matrix.
219:       - Since this is such a small problem, we set all entries for
220:         the matrix at once.
221:   */
222:   A[0]  = 2.0*xx[0] + xx[1]; A[1] = xx[0];
223:   A[2]  = xx[1]; A[3] = xx[0] + 2.0*xx[1];
224:   MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);

226:   /*
227:      Restore vector
228:   */
229:   VecRestoreArrayRead(x,&xx);

231:   /*
232:      Assemble matrix
233:   */
234:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
235:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
236:   if (jac != B) {
237:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
238:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
239:   }
240:   return 0;
241: }

243: /* ------------------------------------------------------------------- */
244: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
245: {
246:   PetscErrorCode    ierr;
247:   const PetscScalar *xx;
248:   PetscScalar       *ff;

250:   /*
251:      Get pointers to vector data.
252:        - For default PETSc vectors, VecGetArray() returns a pointer to
253:          the data array.  Otherwise, the routine is implementation dependent.
254:        - You MUST call VecRestoreArray() when you no longer need access to
255:          the array.
256:   */
257:   VecGetArrayRead(x,&xx);
258:   VecGetArray(f,&ff);

260:   /*
261:      Compute function
262:   */
263:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
264:   ff[1] = xx[1];

266:   /*
267:      Restore vectors
268:   */
269:   VecRestoreArrayRead(x,&xx);
270:   VecRestoreArray(f,&ff);
271:   return 0;
272: }
273: /* ------------------------------------------------------------------- */
274: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
275: {
276:   const PetscScalar *xx;
277:   PetscScalar       A[4];
278:   PetscErrorCode    ierr;
279:   PetscInt          idx[2] = {0,1};

281:   /*
282:      Get pointer to vector data
283:   */
284:   VecGetArrayRead(x,&xx);

286:   /*
287:      Compute Jacobian entries and insert into matrix.
288:       - Since this is such a small problem, we set all entries for
289:         the matrix at once.
290:   */
291:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
292:   A[2]  = 0.0;                     A[3] = 1.0;
293:   MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);

295:   /*
296:      Restore vector
297:   */
298:   VecRestoreArrayRead(x,&xx);

300:   /*
301:      Assemble matrix
302:   */
303:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
304:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
305:   if (jac != B) {
306:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
307:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
308:   }
309:   return 0;
310: }