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*);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

162:    Input Parameters:
163: .  snes - the SNES context
164: .  x    - input vector
165: .  ctx  - optional user-defined context

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

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

186:   /* Compute function */
187:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
188:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

190:   /* Restore vectors */
191:   VecRestoreArrayRead(x,&xx);
192:   VecRestoreArray(f,&ff);
193:   return 0;
194: }
195: /* ------------------------------------------------------------------- */
198: /*
199:    FormJacobian1 - Evaluates Jacobian matrix.

201:    Input Parameters:
202: .  snes - the SNES context
203: .  x - input vector
204: .  dummy - optional user-defined context (not used here)

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

218:   /*
219:      Get pointer to vector data
220:   */
221:   VecGetArrayRead(x,&xx);

223:   /*
224:      Compute Jacobian entries and insert into matrix.
225:       - Since this is such a small problem, we set all entries for
226:         the matrix at once.
227:   */
228:   A[0]  = 2.0*xx[0] + xx[1]; A[1] = xx[0];
229:   A[2]  = xx[1]; A[3] = xx[0] + 2.0*xx[1];
230:   MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);

232:   /*
233:      Restore vector
234:   */
235:   VecRestoreArrayRead(x,&xx);

237:   /*
238:      Assemble matrix
239:   */
240:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
241:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
242:   if (jac != B) {
243:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
244:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
245:   }
246:   return 0;
247: }

249: /* ------------------------------------------------------------------- */
252: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
253: {
254:   PetscErrorCode    ierr;
255:   const PetscScalar *xx;
256:   PetscScalar       *ff;

258:   /*
259:      Get pointers to vector data.
260:        - For default PETSc vectors, VecGetArray() returns a pointer to
261:          the data array.  Otherwise, the routine is implementation dependent.
262:        - You MUST call VecRestoreArray() when you no longer need access to
263:          the array.
264:   */
265:   VecGetArrayRead(x,&xx);
266:   VecGetArray(f,&ff);

268:   /*
269:      Compute function
270:   */
271:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
272:   ff[1] = xx[1];

274:   /*
275:      Restore vectors
276:   */
277:   VecRestoreArrayRead(x,&xx);
278:   VecRestoreArray(f,&ff);
279:   return 0;
280: }
281: /* ------------------------------------------------------------------- */
284: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
285: {
286:   const PetscScalar *xx;
287:   PetscScalar       A[4];
288:   PetscErrorCode    ierr;
289:   PetscInt          idx[2] = {0,1};

291:   /*
292:      Get pointer to vector data
293:   */
294:   VecGetArrayRead(x,&xx);

296:   /*
297:      Compute Jacobian entries and insert into matrix.
298:       - Since this is such a small problem, we set all entries for
299:         the matrix at once.
300:   */
301:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
302:   A[2]  = 0.0;                     A[3] = 1.0;
303:   MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);

305:   /*
306:      Restore vector
307:   */
308:   VecRestoreArrayRead(x,&xx);

310:   /*
311:      Assemble matrix
312:   */
313:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
314:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
315:   if (jac != B) {
316:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
317:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
318:   }
319:   return 0;
320: }