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/examples/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*,MatStructure*,void*);
 49: static PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 50: static PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 51: static PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
 52: static PetscErrorCode FormJacobian1_block(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 53: static PetscErrorCode FormFunction1_block(SNES,Vec,Vec,void*);
 54: static PetscErrorCode FormJacobian2_block(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 55: static PetscErrorCode FormFunction2_block(SNES,Vec,Vec,void*);


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

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

 75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76:   Create nonlinear solver context
 77:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 79:   SNESCreate(PETSC_COMM_WORLD,&snes);

 81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 82:   Create matrix and vector data structures; set corresponding routines
 83:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 85:   /*
 86:   Create vectors for solution and nonlinear function
 87:   */
 88:   VecCreateSeq(PETSC_COMM_SELF,2,&x);
 89:   VecDuplicate(x,&r);

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

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

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

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:   Customize nonlinear solver; set runtime options
116:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

154:   SNESSolve(snes,PETSC_NULL,x);
155:   SNESGetIterationNumber(snes,&its);
156:   if (flg) {
157:     Vec f;
158:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
159:     SNESGetFunction(snes,&f,0,0);
160:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
161:   }

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

165:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166:   Free work space.  All PETSc objects should be destroyed when they
167:   are no longer needed.
168:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

170:   VecDestroy(&x); VecDestroy(&r);
171:   MatDestroy(&J); SNESDestroy(&snes);
172:   return(0);
173: }
174: /* ------------------------------------------------------------------- */
177: /*
178: FormFunction1 - Evaluates nonlinear function, F(x).

180: Input Parameters:
181: .  snes - the SNES context
182: .  x - input vector
183: .  dummy - optional user-defined context (not used here)

185: Output Parameter:
186: .  f - function vector
187: */
188: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
189: {
191:   PetscScalar    *xx,*ff;

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

204:   /*
205:   Compute function
206:   */
207:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
208:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;


211:   /*
212:   Restore vectors
213:   */
214:   VecRestoreArray(x,&xx);
215:   VecRestoreArray(f,&ff);
216:   return(0);
217: }
218: /* ------------------------------------------------------------------- */
221: /*
222: FormJacobian1 - Evaluates Jacobian matrix.

224: Input Parameters:
225: .  snes - the SNES context
226: .  x - input vector
227: .  dummy - optional user-defined context (not used here)

229: Output Parameters:
230: .  jac - Jacobian matrix
231: .  B - optionally different preconditioning matrix
232: .  flag - flag indicating matrix structure
233: */
234: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
235: {
236:   PetscScalar    *xx,A[4];
238:   PetscInt       idx[2] = {0,1};

241:   /*
242:   Get pointer to vector data
243:   */
244:   VecGetArray(x,&xx);

246:   /*
247:   Compute Jacobian entries and insert into matrix.
248:   - Since this is such a small problem, we set all entries for
249:   the matrix at once.
250:   */
251:   A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
252:   A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
253:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
254:   *flag = SAME_NONZERO_PATTERN;

256:   /*
257:   Restore vector
258:   */
259:   VecRestoreArray(x,&xx);

261:   /*
262:   Assemble matrix
263:   */
264:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
265:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
266:   return(0);
267: }


270: /* ------------------------------------------------------------------- */
273: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
274: {
276:   PetscScalar    *xx,*ff;

279:   /*
280:   Get pointers to vector data.
281:   - For default PETSc vectors, VecGetArray() returns a pointer to
282:   the data array.  Otherwise, the routine is implementation dependent.
283:   - You MUST call VecRestoreArray() when you no longer need access to
284:   the array.
285:   */
286:   VecGetArray(x,&xx);
287:   VecGetArray(f,&ff);

289:   /*
290:   Compute function
291:   */
292:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
293:   ff[1] = xx[1];

295:   /*
296:   Restore vectors
297:   */
298:   VecRestoreArray(x,&xx);
299:   VecRestoreArray(f,&ff);
300:   return(0);
301: }
302: /* ------------------------------------------------------------------- */
305: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
306: {
307:   PetscScalar    *xx,A[4];
309:   PetscInt       idx[2] = {0,1};

312:   /*
313:   Get pointer to vector data
314:   */
315:   VecGetArray(x,&xx);

317:   /*
318:   Compute Jacobian entries and insert into matrix.
319:   - Since this is such a small problem, we set all entries for
320:   the matrix at once.
321:   */
322:   A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
323:   A[2] = 0.0;                     A[3] = 1.0;
324:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
325:   *flag = SAME_NONZERO_PATTERN;

327:   /*
328:   Restore vector
329:   */
330:   VecRestoreArray(x,&xx);

332:   /*
333:   Assemble matrix
334:   */
335:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
336:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
337:   return(0);
338: }

342: static int block_system(void)
343: {
344:   SNES           snes;         /* nonlinear solver context */
345:   KSP            ksp;         /* linear solver context */
346:   PC             pc;           /* preconditioner context */
347:   Vec            x,r;         /* solution, residual vectors */
348:   Mat            J;            /* Jacobian matrix */
350:   PetscInt       its;
351:   PetscScalar    pfive = .5;
352:   PetscBool      flg;

354:   Mat j11, j12, j21, j22;
355:   Vec x1, x2, r1, r2;
356:   Vec bv;
357:   Vec bx[2];
358:   Mat bA[2][2];

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

363:   SNESCreate(PETSC_COMM_WORLD,&snes);

365:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
366:   Create matrix and vector data structures; set corresponding routines
367:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

369:   /*
370:   Create sub vectors for solution and nonlinear function
371:   */
372:   VecCreateSeq(PETSC_COMM_SELF,1,&x1);
373:   VecDuplicate(x1,&r1);

375:   VecCreateSeq(PETSC_COMM_SELF,1,&x2);
376:   VecDuplicate(x2,&r2);

378:   /*
379:   Create the block vectors
380:   */
381:   bx[0] = x1;
382:   bx[1] = x2;
383:   VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,bx,&x);
384:   VecAssemblyBegin(x);
385:   VecAssemblyEnd(x);
386:   VecDestroy(&x1);
387:   VecDestroy(&x2);

389:   bx[0] = r1;
390:   bx[1] = r2;
391:   VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,bx,&r);
392:   VecDestroy(&r1);
393:   VecDestroy(&r2);
394:   VecAssemblyBegin(r);
395:   VecAssemblyEnd(r);

397:   /*
398:   Create sub Jacobian matrix data structure
399:   */
400:   MatCreate( PETSC_COMM_WORLD, &j11 );
401:   MatSetSizes( j11, 1, 1, 1, 1 );
402:   MatSetType( j11, MATSEQAIJ );

404:   MatCreate( PETSC_COMM_WORLD, &j12 );
405:   MatSetSizes( j12, 1, 1, 1, 1 );
406:   MatSetType( j12, MATSEQAIJ );

408:   MatCreate( PETSC_COMM_WORLD, &j21 );
409:   MatSetSizes( j21, 1, 1, 1, 1 );
410:   MatSetType( j21, MATSEQAIJ );

412:   MatCreate( PETSC_COMM_WORLD, &j22 );
413:   MatSetSizes( j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1 );
414:   MatSetType( j22, MATSEQAIJ );

416:   /*
417:   Create block Jacobian matrix data structure
418:   */
419:   bA[0][0] = j11;
420:   bA[0][1] = j12;
421:   bA[1][0] = j21;
422:   bA[1][1] = j22;
423:   MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,2,PETSC_NULL,&bA[0][0],&J);
424:   MatNestSetVecType(J,VECNEST);
425:   MatDestroy(&j11);
426:   MatDestroy(&j12);
427:   MatDestroy(&j21);
428:   MatDestroy(&j22);

430:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
431:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

433:   PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
434:   if (!flg) {
435:     /*
436:     Set function evaluation routine and vector.
437:     */
438:     SNESSetFunction(snes,r,FormFunction1_block,PETSC_NULL);

440:     /*
441:     Set Jacobian matrix data structure and Jacobian evaluation routine
442:     */
443:     SNESSetJacobian(snes,J,J,FormJacobian1_block,PETSC_NULL);
444:   } else {
445:     SNESSetFunction(snes,r,FormFunction2_block,PETSC_NULL);
446:     SNESSetJacobian(snes,J,J,FormJacobian2_block,PETSC_NULL);
447:   }

449:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
450:   Customize nonlinear solver; set runtime options
451:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

463:   /*
464:   Set SNES/KSP/KSP/PC runtime options, e.g.,
465:   -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
466:   These options will override those specified above as long as
467:   SNESSetFromOptions() is called _after_ any other customization
468:   routines.
469:   */
470:   SNESSetFromOptions(snes);

472:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
473:   Evaluate initial guess; then solve nonlinear system
474:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
475:   if (!flg) {
476:     VecSet(x,pfive);
477:   } else {
478:     Vec *vecs;
479:     VecNestGetSubVecs( x, PETSC_NULL, &vecs );
480:     bv = vecs[0];
481: //    VecBlockGetSubVec( x, 0, &bv );
482:     VecSetValue( bv, 0, 2.0, INSERT_VALUES );  /* xx[0] = 2.0; */
483:     VecAssemblyBegin(bv);
484:     VecAssemblyEnd(bv);

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

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

509:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
510:   Free work space.  All PETSc objects should be destroyed when they
511:   are no longer needed.
512:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
513:   VecDestroy(&x); VecDestroy(&r);
514:   MatDestroy(&J); SNESDestroy(&snes);

516:   return(0);
517: }
518: /* ------------------------------------------------------------------- */
521: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
522: {
523:   Vec *xx, *ff, x1,x2, f1,f2;
524:   PetscScalar ff_0, ff_1;
525:   PetscScalar xx_0, xx_1;
526:   PetscInt index,nb;

530:   /* get blocks for function */
531:   VecNestGetSubVecs( f, &nb, &ff );
532:   f1 = ff[0];  f2 = ff[1];

534:   /* get blocks for solution */
535:   VecNestGetSubVecs( x, &nb, &xx );
536:   x1 = xx[0];  x2 = xx[1];

538:   /* get solution values */
539:   index = 0;
540:   VecGetValues( x1,1, &index, &xx_0 );
541:   VecGetValues( x2,1, &index, &xx_1 );

543:   /* Compute function */
544:   ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
545:   ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;

547:   /* set function values */
548:   VecSetValue( f1, index, ff_0, INSERT_VALUES );

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

552:   VecAssemblyBegin(f);
553:   VecAssemblyEnd(f);

555:   return(0);
556: }
557: /* ------------------------------------------------------------------- */
560: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
561: {
562:   Vec            *xx, x1,x2;
563:   PetscScalar    xx_0, xx_1;
564:   PetscInt       index,nb;
565:   PetscScalar    A_00, A_01, A_10, A_11;
566:   Mat            j11, j12, j21, j22;
567:   Mat            **mats;

571:   /* get blocks for solution */
572:   VecNestGetSubVecs( x, &nb, &xx );
573:   x1 = xx[0];  x2 = xx[1];

575:   /* get solution values */
576:   index = 0;
577:   VecGetValues( x1,1, &index, &xx_0 );
578:   VecGetValues( x2,1, &index, &xx_1 );

580:   /* get block matrices */
581:   MatNestGetSubMats(*jac,PETSC_NULL,PETSC_NULL,&mats);
582:   j11 = mats[0][0];
583:   j12 = mats[0][1];
584:   j21 = mats[1][0];
585:   j22 = mats[1][1];

587:   /* compute jacobian entries */
588:   A_00 = 2.0*xx_0 + xx_1;
589:   A_01 = xx_0;
590:   A_10 = xx_1;
591:   A_11 = xx_0 + 2.0*xx_1;

593:   /* set jacobian values */
594:   MatSetValue( j11, 0,0, A_00, INSERT_VALUES);
595:   MatSetValue( j12, 0,0, A_01, INSERT_VALUES);
596:   MatSetValue( j21, 0,0, A_10, INSERT_VALUES);
597:   MatSetValue( j22, 0,0, A_11, INSERT_VALUES);

599:   *flag = SAME_NONZERO_PATTERN;

601:   /* Assemble sub matrix */
602:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
603:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
604:   return(0);
605: }


608: /* ------------------------------------------------------------------- */
611: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
612: {
614:   PetscScalar    *xx,*ff;

617:   /*
618:   Get pointers to vector data.
619:   - For default PETSc vectors, VecGetArray() returns a pointer to
620:   the data array.  Otherwise, the routine is implementation dependent.
621:   - You MUST call VecRestoreArray() when you no longer need access to
622:   the array.
623:   */
624:   VecGetArray(x,&xx);
625:   VecGetArray(f,&ff);

627:   /*
628:   Compute function
629:   */
630:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
631:   ff[1] = xx[1];

633:   /*
634:   Restore vectors
635:   */
636:   VecRestoreArray(x,&xx);
637:   VecRestoreArray(f,&ff);
638:   return(0);
639: }
640: /* ------------------------------------------------------------------- */
643: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
644: {
645:   PetscScalar    *xx,A[4];
647:   PetscInt       idx[2] = {0,1};

650:   /*
651:   Get pointer to vector data
652:   */
653:   VecGetArray(x,&xx);

655:   /*
656:   Compute Jacobian entries and insert into matrix.
657:   - Since this is such a small problem, we set all entries for
658:   the matrix at once.
659:   */
660:   A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
661:   A[2] = 0.0;                     A[3] = 1.0;
662:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
663:   *flag = SAME_NONZERO_PATTERN;

665:   /*
666:   Restore vector
667:   */
668:   VecRestoreArray(x,&xx);

670:   /*
671:   Assemble matrix
672:   */
673:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
674:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
675:   return(0);
676: }

680: int main(int argc,char **argv)
681: {
682:   PetscMPIInt    size;

685:   PetscInitialize(&argc,&argv,(char *)0,help);

687:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
688:   if (size != 1) SETERRQ(PETSC_COMM_WORLD, 1,"This is a uniprocessor example only!");

690:   assembled_system();

692:   block_system();

694:   PetscFinalize();
695:   return 0;
696: }