Actual source code: ex17.c

  1: static const char help[] = "Newton's method to solve a two-variable system, sequentially.\n"
  2:                            "The same problem is solved twice - i) fully assembled system + ii) block system\n\n";

  4: /*T
  5: Concepts: SNES^basic uniprocessor example, block objects
  6: Processors: 1
  7: T*/

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

 20: /*
 21: This example is block version of the test found at
 22:   ${PETSC_DIR}/src/snes/tutorials/ex1.c
 23: In this test we replace the Jacobian systems
 24:   [J]{x} = {F}
 25: where

 27: [J] = (j_00, j_01),  {x} = (x_0, x_1)^T,   {F} = (f_0, f_1)^T
 28:       (j_10, j_11)
 29: where [J] \in \mathbb^{2 \times 2}, {x},{F} \in \mathbb^{2 \times 1},

 31: with a block system in which each block is of length 1.
 32: i.e. The block system is thus

 34: [J] = ([j00], [j01]),  {x} = ({x0}, {x1})^T, {F} = ({f0}, {f1})^T
 35:       ([j10], [j11])
 36: where
 37: [j00], [j01], [j10], [j11] \in \mathbb^{1 \times 1}
 38: {x0}, {x1}, {f0}, {f1} \in \mathbb^{1 \times 1}

 40: In practice we would not bother defing blocks of size one, and would instead assemble the
 41: full system. This is just a simple test to illustrate how to manipulate the blocks and
 42: to confirm the implementation is correct.
 43: */

 45: /*
 46: User-defined routines
 47: */
 48: static PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
 49: static PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 50: static PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
 51: static PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
 52: static PetscErrorCode FormJacobian1_block(SNES,Vec,Mat,Mat,void*);
 53: static PetscErrorCode FormFunction1_block(SNES,Vec,Vec,void*);
 54: static PetscErrorCode FormJacobian2_block(SNES,Vec,Mat,Mat,void*);
 55: static PetscErrorCode FormFunction2_block(SNES,Vec,Vec,void*);

 57: static PetscErrorCode assembled_system(void)
 58: {
 59:   SNES           snes;         /* nonlinear solver context */
 60:   KSP            ksp;         /* linear solver context */
 61:   PC             pc;           /* preconditioner context */
 62:   Vec            x,r;         /* solution, residual vectors */
 63:   Mat            J;            /* Jacobian matrix */
 65:   PetscInt       its;
 66:   PetscScalar    pfive = .5,*xx;
 67:   PetscBool      flg;

 70:   PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n");

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 73:   Create nonlinear solver context
 74:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 76:   SNESCreate(PETSC_COMM_WORLD,&snes);

 78:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79:   Create matrix and vector data structures; set corresponding routines
 80:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 82:   /*
 83:   Create vectors for solution and nonlinear function
 84:   */
 85:   VecCreateSeq(PETSC_COMM_SELF,2,&x);
 86:   VecDuplicate(x,&r);

 88:   /*
 89:   Create Jacobian matrix data structure
 90:   */
 91:   MatCreate(PETSC_COMM_SELF,&J);
 92:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 93:   MatSetFromOptions(J);
 94:   MatSetUp(J);

 96:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
 97:   if (!flg) {
 98:     /*
 99:     Set function evaluation routine and vector.
100:     */
101:     SNESSetFunction(snes,r,FormFunction1,NULL);

103:     /*
104:     Set Jacobian matrix data structure and Jacobian evaluation routine
105:     */
106:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
107:   } else {
108:     SNESSetFunction(snes,r,FormFunction2,NULL);
109:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
110:   }

112:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113:   Customize nonlinear solver; set runtime options
114:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

116:   /*
117:   Set linear solver defaults for this problem. By extracting the
118:   KSP, KSP, and PC contexts from the SNES context, we can then
119:   directly call any KSP, KSP, and PC routines to set various options.
120:   */
121:   SNESGetKSP(snes,&ksp);
122:   KSPGetPC(ksp,&pc);
123:   PCSetType(pc,PCNONE);
124:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

126:   /*
127:   Set SNES/KSP/KSP/PC runtime options, e.g.,
128:   -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
129:   These options will override those specified above as long as
130:   SNESSetFromOptions() is called _after_ any other customization
131:   routines.
132:   */
133:   SNESSetFromOptions(snes);

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:   Evaluate initial guess; then solve nonlinear system
137:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138:   if (!flg) {
139:     VecSet(x,pfive);
140:   } else {
141:     VecGetArray(x,&xx);
142:     xx[0] = 2.0; xx[1] = 3.0;
143:     VecRestoreArray(x,&xx);
144:   }
145:   /*
146:   Note: The user should initialize the vector, x, with the initial guess
147:   for the nonlinear solver prior to calling SNESSolve().  In particular,
148:   to employ an initial guess of zero, the user should explicitly set
149:   this vector to zero by calling VecSet().
150:   */

152:   SNESSolve(snes,NULL,x);
153:   SNESGetIterationNumber(snes,&its);
154:   if (flg) {
155:     Vec f;
156:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
157:     SNESGetFunction(snes,&f,0,0);
158:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
159:   }
160:   PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);

162:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163:   Free work space.  All PETSc objects should be destroyed when they
164:   are no longer needed.
165:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

167:   VecDestroy(&x); VecDestroy(&r);
168:   MatDestroy(&J); SNESDestroy(&snes);
169:   return(0);
170: }

172: /* ------------------------------------------------------------------- */
173: /*
174: FormFunction1 - Evaluates nonlinear function, F(x).

176: Input Parameters:
177: .  snes - the SNES context
178: .  x - input vector
179: .  dummy - optional user-defined context (not used here)

181: Output Parameter:
182: .  f - function vector
183: */
184: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
185: {
186:   PetscErrorCode    ierr;
187:   const PetscScalar *xx;
188:   PetscScalar       *ff;

191:   /*
192:   Get pointers to vector data.
193:   - For default PETSc vectors, VecGetArray() returns a pointer to
194:   the data array.  Otherwise, the routine is implementation dependent.
195:   - You MUST call VecRestoreArray() when you no longer need access to
196:   the array.
197:   */
198:   VecGetArrayRead(x,&xx);
199:   VecGetArray(f,&ff);

201:   /*
202:   Compute function
203:   */
204:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
205:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

207:   /*
208:   Restore vectors
209:   */
210:   VecRestoreArrayRead(x,&xx);
211:   VecRestoreArray(f,&ff);
212:   return(0);
213: }
214: /* ------------------------------------------------------------------- */
215: /*
216: FormJacobian1 - Evaluates Jacobian matrix.

218: Input Parameters:
219: .  snes - the SNES context
220: .  x - input vector
221: .  dummy - optional user-defined context (not used here)

223: Output Parameters:
224: .  jac - Jacobian matrix
225: .  B - optionally different preconditioning matrix
226: .  flag - flag indicating matrix structure
227: */
228: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
229: {
230:   const PetscScalar *xx;
231:   PetscScalar       A[4];
232:   PetscErrorCode    ierr;
233:   PetscInt          idx[2] = {0,1};

236:   /*
237:   Get pointer to vector data
238:   */
239:   VecGetArrayRead(x,&xx);

241:   /*
242:   Compute Jacobian entries and insert into matrix.
243:   - Since this is such a small problem, we set all entries for
244:   the matrix at once.
245:   */
246:   A[0]  = 2.0*xx[0] + xx[1]; A[1] = xx[0];
247:   A[2]  = xx[1]; A[3] = xx[0] + 2.0*xx[1];
248:   MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);

250:   /*
251:   Restore vector
252:   */
253:   VecRestoreArrayRead(x,&xx);

255:   /*
256:   Assemble matrix
257:   */
258:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
259:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
260:   return(0);
261: }

263: /* ------------------------------------------------------------------- */
264: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
265: {
266:   PetscErrorCode    ierr;
267:   const PetscScalar *xx;
268:   PetscScalar       *ff;

271:   /*
272:   Get pointers to vector data.
273:   - For default PETSc vectors, VecGetArray() returns a pointer to
274:   the data array.  Otherwise, the routine is implementation dependent.
275:   - You MUST call VecRestoreArray() when you no longer need access to
276:   the array.
277:   */
278:   VecGetArrayRead(x,&xx);
279:   VecGetArray(f,&ff);

281:   /*
282:   Compute function
283:   */
284:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
285:   ff[1] = xx[1];

287:   /*
288:   Restore vectors
289:   */
290:   VecRestoreArrayRead(x,&xx);
291:   VecRestoreArray(f,&ff);
292:   return(0);
293: }

295: /* ------------------------------------------------------------------- */
296: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
297: {
298:   const PetscScalar *xx;
299:   PetscScalar       A[4];
300:   PetscErrorCode    ierr;
301:   PetscInt          idx[2] = {0,1};

304:   /*
305:   Get pointer to vector data
306:   */
307:   VecGetArrayRead(x,&xx);

309:   /*
310:   Compute Jacobian entries and insert into matrix.
311:   - Since this is such a small problem, we set all entries for
312:   the matrix at once.
313:   */
314:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
315:   A[2]  = 0.0;                     A[3] = 1.0;
316:   MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);

318:   /*
319:   Restore vector
320:   */
321:   VecRestoreArrayRead(x,&xx);

323:   /*
324:   Assemble matrix
325:   */
326:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
327:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
328:   return(0);
329: }

331: static PetscErrorCode block_system(void)
332: {
333:   SNES           snes;         /* nonlinear solver context */
334:   KSP            ksp;         /* linear solver context */
335:   PC             pc;           /* preconditioner context */
336:   Vec            x,r;         /* solution, residual vectors */
337:   Mat            J;            /* Jacobian matrix */
339:   PetscInt       its;
340:   PetscScalar    pfive = .5;
341:   PetscBool      flg;

343:   Mat            j11, j12, j21, j22;
344:   Vec            x1, x2, r1, r2;
345:   Vec            bv;
346:   Vec            bx[2];
347:   Mat            bA[2][2];

350:   PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n");

352:   SNESCreate(PETSC_COMM_WORLD,&snes);

354:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
355:   Create matrix and vector data structures; set corresponding routines
356:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

358:   /*
359:   Create sub vectors for solution and nonlinear function
360:   */
361:   VecCreateSeq(PETSC_COMM_SELF,1,&x1);
362:   VecDuplicate(x1,&r1);

364:   VecCreateSeq(PETSC_COMM_SELF,1,&x2);
365:   VecDuplicate(x2,&r2);

367:   /*
368:   Create the block vectors
369:   */
370:   bx[0] = x1;
371:   bx[1] = x2;
372:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
373:   VecAssemblyBegin(x);
374:   VecAssemblyEnd(x);
375:   VecDestroy(&x1);
376:   VecDestroy(&x2);

378:   bx[0] = r1;
379:   bx[1] = r2;
380:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
381:   VecDestroy(&r1);
382:   VecDestroy(&r2);
383:   VecAssemblyBegin(r);
384:   VecAssemblyEnd(r);

386:   /*
387:   Create sub Jacobian matrix data structure
388:   */
389:   MatCreate(PETSC_COMM_WORLD, &j11);
390:   MatSetSizes(j11, 1, 1, 1, 1);
391:   MatSetType(j11, MATSEQAIJ);
392:   MatSetUp(j11);

394:   MatCreate(PETSC_COMM_WORLD, &j12);
395:   MatSetSizes(j12, 1, 1, 1, 1);
396:   MatSetType(j12, MATSEQAIJ);
397:   MatSetUp(j12);

399:   MatCreate(PETSC_COMM_WORLD, &j21);
400:   MatSetSizes(j21, 1, 1, 1, 1);
401:   MatSetType(j21, MATSEQAIJ);
402:   MatSetUp(j21);

404:   MatCreate(PETSC_COMM_WORLD, &j22);
405:   MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
406:   MatSetType(j22, MATSEQAIJ);
407:   MatSetUp(j22);
408:   /*
409:   Create block Jacobian matrix data structure
410:   */
411:   bA[0][0] = j11;
412:   bA[0][1] = j12;
413:   bA[1][0] = j21;
414:   bA[1][1] = j22;

416:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
417:   MatSetUp(J);
418:   MatNestSetVecType(J,VECNEST);
419:   MatDestroy(&j11);
420:   MatDestroy(&j12);
421:   MatDestroy(&j21);
422:   MatDestroy(&j22);

424:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
425:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

427:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
428:   if (!flg) {
429:     /*
430:     Set function evaluation routine and vector.
431:     */
432:     SNESSetFunction(snes,r,FormFunction1_block,NULL);

434:     /*
435:     Set Jacobian matrix data structure and Jacobian evaluation routine
436:     */
437:     SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
438:   } else {
439:     SNESSetFunction(snes,r,FormFunction2_block,NULL);
440:     SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
441:   }

443:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
444:   Customize nonlinear solver; set runtime options
445:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

447:   /*
448:   Set linear solver defaults for this problem. By extracting the
449:   KSP, KSP, and PC contexts from the SNES context, we can then
450:   directly call any KSP, KSP, and PC routines to set various options.
451:   */
452:   SNESGetKSP(snes,&ksp);
453:   KSPGetPC(ksp,&pc);
454:   PCSetType(pc,PCNONE);
455:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

457:   /*
458:   Set SNES/KSP/KSP/PC runtime options, e.g.,
459:   -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
460:   These options will override those specified above as long as
461:   SNESSetFromOptions() is called _after_ any other customization
462:   routines.
463:   */
464:   SNESSetFromOptions(snes);

466:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
467:   Evaluate initial guess; then solve nonlinear system
468:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
469:   if (!flg) {
470:     VecSet(x,pfive);
471:   } else {
472:     Vec *vecs;
473:     VecNestGetSubVecs(x, NULL, &vecs);
474:     bv   = vecs[0];
475: /*    VecBlockGetSubVec(x, 0, &bv); */
476:     VecSetValue(bv, 0, 2.0, INSERT_VALUES);  /* xx[0] = 2.0; */
477:     VecAssemblyBegin(bv);
478:     VecAssemblyEnd(bv);

480: /*    VecBlockGetSubVec(x, 1, &bv); */
481:     bv   = vecs[1];
482:     VecSetValue(bv, 0, 3.0, INSERT_VALUES);  /* xx[1] = 3.0; */
483:     VecAssemblyBegin(bv);
484:     VecAssemblyEnd(bv);
485:   }
486:   /*
487:   Note: The user should initialize the vector, x, with the initial guess
488:   for the nonlinear solver prior to calling SNESSolve().  In particular,
489:   to employ an initial guess of zero, the user should explicitly set
490:   this vector to zero by calling VecSet().
491:   */
492:   SNESSolve(snes,NULL,x);
493:   SNESGetIterationNumber(snes,&its);
494:   if (flg) {
495:     Vec f;
496:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
497:     SNESGetFunction(snes,&f,0,0);
498:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
499:   }

501:   PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);

503:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
504:   Free work space.  All PETSc objects should be destroyed when they
505:   are no longer needed.
506:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
507:   VecDestroy(&x); VecDestroy(&r);
508:   MatDestroy(&J); SNESDestroy(&snes);
509:   return(0);
510: }

512: /* ------------------------------------------------------------------- */
513: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
514: {
515:   Vec            *xx, *ff, x1,x2, f1,f2;
516:   PetscScalar    ff_0, ff_1;
517:   PetscScalar    xx_0, xx_1;
518:   PetscInt       index,nb;

522:   /* get blocks for function */
523:   VecNestGetSubVecs(f, &nb, &ff);
524:   f1   = ff[0];  f2 = ff[1];

526:   /* get blocks for solution */
527:   VecNestGetSubVecs(x, &nb, &xx);
528:   x1   = xx[0];  x2 = xx[1];

530:   /* get solution values */
531:   index = 0;
532:   VecGetValues(x1,1, &index, &xx_0);
533:   VecGetValues(x2,1, &index, &xx_1);

535:   /* Compute function */
536:   ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
537:   ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;

539:   /* set function values */
540:   VecSetValue(f1, index, ff_0, INSERT_VALUES);

542:   VecSetValue(f2, index, ff_1, INSERT_VALUES);

544:   VecAssemblyBegin(f);
545:   VecAssemblyEnd(f);
546:   return(0);
547: }

549: /* ------------------------------------------------------------------- */
550: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
551: {
552:   Vec            *xx, x1,x2;
553:   PetscScalar    xx_0, xx_1;
554:   PetscInt       index,nb;
555:   PetscScalar    A_00, A_01, A_10, A_11;
556:   Mat            j11, j12, j21, j22;
557:   Mat            **mats;

561:   /* get blocks for solution */
562:   VecNestGetSubVecs(x, &nb, &xx);
563:   x1   = xx[0];  x2 = xx[1];

565:   /* get solution values */
566:   index = 0;
567:   VecGetValues(x1,1, &index, &xx_0);
568:   VecGetValues(x2,1, &index, &xx_1);

570:   /* get block matrices */
571:   MatNestGetSubMats(jac,NULL,NULL,&mats);
572:   j11  = mats[0][0];
573:   j12  = mats[0][1];
574:   j21  = mats[1][0];
575:   j22  = mats[1][1];

577:   /* compute jacobian entries */
578:   A_00 = 2.0*xx_0 + xx_1;
579:   A_01 = xx_0;
580:   A_10 = xx_1;
581:   A_11 = xx_0 + 2.0*xx_1;

583:   /* set jacobian values */
584:   MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
585:   MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
586:   MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
587:   MatSetValue(j22, 0,0, A_11, INSERT_VALUES);

589:   /* Assemble sub matrix */
590:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
591:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
592:   return(0);
593: }

595: /* ------------------------------------------------------------------- */
596: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
597: {
598:   PetscErrorCode    ierr;
599:   PetscScalar       *ff;
600:   const PetscScalar *xx;

603:   /*
604:   Get pointers to vector data.
605:   - For default PETSc vectors, VecGetArray() returns a pointer to
606:   the data array.  Otherwise, the routine is implementation dependent.
607:   - You MUST call VecRestoreArray() when you no longer need access to
608:   the array.
609:   */
610:   VecGetArrayRead(x,&xx);
611:   VecGetArray(f,&ff);

613:   /*
614:   Compute function
615:   */
616:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
617:   ff[1] = xx[1];

619:   /*
620:   Restore vectors
621:   */
622:   VecRestoreArrayRead(x,&xx);
623:   VecRestoreArray(f,&ff);
624:   return(0);
625: }

627: /* ------------------------------------------------------------------- */
628: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
629: {
630:   const PetscScalar *xx;
631:   PetscScalar       A[4];
632:   PetscErrorCode    ierr;
633:   PetscInt          idx[2] = {0,1};

636:   /*
637:   Get pointer to vector data
638:   */
639:   VecGetArrayRead(x,&xx);

641:   /*
642:   Compute Jacobian entries and insert into matrix.
643:   - Since this is such a small problem, we set all entries for
644:   the matrix at once.
645:   */
646:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
647:   A[2]  = 0.0;                     A[3] = 1.0;
648:   MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);

650:   /*
651:   Restore vector
652:   */
653:   VecRestoreArrayRead(x,&xx);

655:   /*
656:   Assemble matrix
657:   */
658:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
659:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
660:   return(0);
661: }

663: int main(int argc,char **argv)
664: {
665:   PetscMPIInt    size;

668:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
669:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
670:   if (size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
671:   assembled_system();
672:   block_system();
673:   PetscFinalize();
674:   return ierr;
675: }

677: /*TEST

679:    test:
680:       args: -snes_monitor_short
681:       requires: !single

683: TEST*/