Actual source code: ex17.c

petsc-3.3-p7 2013-05-11
  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);
 97:   MatSetUp(J);

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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


271: /* ------------------------------------------------------------------- */
274: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
275: {
277:   PetscScalar    *xx,*ff;

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

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

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

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

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

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

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

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

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

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

364:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

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

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

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

406:   MatCreate( PETSC_COMM_WORLD, &j12 );
407:   MatSetSizes( j12, 1, 1, 1, 1 );
408:   MatSetType( j12, MATSEQAIJ );
409:   MatSetUp(j12);

411:   MatCreate( PETSC_COMM_WORLD, &j21 );
412:   MatSetSizes( j21, 1, 1, 1, 1 );
413:   MatSetType( j21, MATSEQAIJ );
414:   MatSetUp(j21);

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

435:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
436:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

438:   PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
439:   if (!flg) {
440:     /*
441:     Set function evaluation routine and vector.
442:     */
443:     SNESSetFunction(snes,r,FormFunction1_block,PETSC_NULL);

445:     /*
446:     Set Jacobian matrix data structure and Jacobian evaluation routine
447:     */
448:     SNESSetJacobian(snes,J,J,FormJacobian1_block,PETSC_NULL);
449:   } else {
450:     SNESSetFunction(snes,r,FormFunction2_block,PETSC_NULL);
451:     SNESSetJacobian(snes,J,J,FormJacobian2_block,PETSC_NULL);
452:   }

454:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
455:   Customize nonlinear solver; set runtime options
456:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

514:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
515:   Free work space.  All PETSc objects should be destroyed when they
516:   are no longer needed.
517:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
518:   VecDestroy(&x); VecDestroy(&r);
519:   MatDestroy(&J); SNESDestroy(&snes);

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

535:   /* get blocks for function */
536:   VecNestGetSubVecs( f, &nb, &ff );
537:   f1 = ff[0];  f2 = ff[1];

539:   /* get blocks for solution */
540:   VecNestGetSubVecs( x, &nb, &xx );
541:   x1 = xx[0];  x2 = xx[1];

543:   /* get solution values */
544:   index = 0;
545:   VecGetValues( x1,1, &index, &xx_0 );
546:   VecGetValues( x2,1, &index, &xx_1 );

548:   /* Compute function */
549:   ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
550:   ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;

552:   /* set function values */
553:   VecSetValue( f1, index, ff_0, INSERT_VALUES );

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

557:   VecAssemblyBegin(f);
558:   VecAssemblyEnd(f);

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

576:   /* get blocks for solution */
577:   VecNestGetSubVecs( x, &nb, &xx );
578:   x1 = xx[0];  x2 = xx[1];

580:   /* get solution values */
581:   index = 0;
582:   VecGetValues( x1,1, &index, &xx_0 );
583:   VecGetValues( x2,1, &index, &xx_1 );

585:   /* get block matrices */
586:   MatNestGetSubMats(*jac,PETSC_NULL,PETSC_NULL,&mats);
587:   j11 = mats[0][0];
588:   j12 = mats[0][1];
589:   j21 = mats[1][0];
590:   j22 = mats[1][1];

592:   /* compute jacobian entries */
593:   A_00 = 2.0*xx_0 + xx_1;
594:   A_01 = xx_0;
595:   A_10 = xx_1;
596:   A_11 = xx_0 + 2.0*xx_1;

598:   /* set jacobian values */
599:   MatSetValue( j11, 0,0, A_00, INSERT_VALUES);
600:   MatSetValue( j12, 0,0, A_01, INSERT_VALUES);
601:   MatSetValue( j21, 0,0, A_10, INSERT_VALUES);
602:   MatSetValue( j22, 0,0, A_11, INSERT_VALUES);

604:   *flag = SAME_NONZERO_PATTERN;

606:   /* Assemble sub matrix */
607:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
608:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
609:   return(0);
610: }


613: /* ------------------------------------------------------------------- */
616: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
617: {
619:   PetscScalar    *xx,*ff;

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

632:   /*
633:   Compute function
634:   */
635:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
636:   ff[1] = xx[1];

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

655:   /*
656:   Get pointer to vector data
657:   */
658:   VecGetArray(x,&xx);

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

670:   /*
671:   Restore vector
672:   */
673:   VecRestoreArray(x,&xx);

675:   /*
676:   Assemble matrix
677:   */
678:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
679:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
680:   return(0);
681: }

685: int main(int argc,char **argv)
686: {
687:   PetscMPIInt    size;

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

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

695:   assembled_system();

697:   block_system();

699:   PetscFinalize();
700:   return 0;
701: }