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 */
 64:   PetscInt       its;
 65:   PetscScalar    pfive = .5,*xx;
 66:   PetscBool      flg;

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

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

 75:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

233:   /*
234:   Get pointer to vector data
235:   */
236:   VecGetArrayRead(x,&xx);

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

247:   /*
248:   Restore vector
249:   */
250:   VecRestoreArrayRead(x,&xx);

252:   /*
253:   Assemble matrix
254:   */
255:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
256:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
257:   return 0;
258: }

260: /* ------------------------------------------------------------------- */
261: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
262: {
263:   const PetscScalar *xx;
264:   PetscScalar       *ff;

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

277:   /*
278:   Compute function
279:   */
280:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
281:   ff[1] = xx[1];

283:   /*
284:   Restore vectors
285:   */
286:   VecRestoreArrayRead(x,&xx);
287:   VecRestoreArray(f,&ff);
288:   return 0;
289: }

291: /* ------------------------------------------------------------------- */
292: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
293: {
294:   const PetscScalar *xx;
295:   PetscScalar       A[4];
296:   PetscInt          idx[2] = {0,1};

299:   /*
300:   Get pointer to vector data
301:   */
302:   VecGetArrayRead(x,&xx);

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

313:   /*
314:   Restore vector
315:   */
316:   VecRestoreArrayRead(x,&xx);

318:   /*
319:   Assemble matrix
320:   */
321:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
322:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
323:   return 0;
324: }

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

337:   Mat            j11, j12, j21, j22;
338:   Vec            x1, x2, r1, r2;
339:   Vec            bv;
340:   Vec            bx[2];
341:   Mat            bA[2][2];

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

346:   SNESCreate(PETSC_COMM_WORLD,&snes);

348:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349:   Create matrix and vector data structures; set corresponding routines
350:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

352:   /*
353:   Create sub vectors for solution and nonlinear function
354:   */
355:   VecCreateSeq(PETSC_COMM_SELF,1,&x1);
356:   VecDuplicate(x1,&r1);

358:   VecCreateSeq(PETSC_COMM_SELF,1,&x2);
359:   VecDuplicate(x2,&r2);

361:   /*
362:   Create the block vectors
363:   */
364:   bx[0] = x1;
365:   bx[1] = x2;
366:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
367:   VecAssemblyBegin(x);
368:   VecAssemblyEnd(x);
369:   VecDestroy(&x1);
370:   VecDestroy(&x2);

372:   bx[0] = r1;
373:   bx[1] = r2;
374:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
375:   VecDestroy(&r1);
376:   VecDestroy(&r2);
377:   VecAssemblyBegin(r);
378:   VecAssemblyEnd(r);

380:   /*
381:   Create sub Jacobian matrix data structure
382:   */
383:   MatCreate(PETSC_COMM_WORLD, &j11);
384:   MatSetSizes(j11, 1, 1, 1, 1);
385:   MatSetType(j11, MATSEQAIJ);
386:   MatSetUp(j11);

388:   MatCreate(PETSC_COMM_WORLD, &j12);
389:   MatSetSizes(j12, 1, 1, 1, 1);
390:   MatSetType(j12, MATSEQAIJ);
391:   MatSetUp(j12);

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

398:   MatCreate(PETSC_COMM_WORLD, &j22);
399:   MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
400:   MatSetType(j22, MATSEQAIJ);
401:   MatSetUp(j22);
402:   /*
403:   Create block Jacobian matrix data structure
404:   */
405:   bA[0][0] = j11;
406:   bA[0][1] = j12;
407:   bA[1][0] = j21;
408:   bA[1][1] = j22;

410:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
411:   MatSetUp(J);
412:   MatNestSetVecType(J,VECNEST);
413:   MatDestroy(&j11);
414:   MatDestroy(&j12);
415:   MatDestroy(&j21);
416:   MatDestroy(&j22);

418:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
419:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

421:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
422:   if (!flg) {
423:     /*
424:     Set function evaluation routine and vector.
425:     */
426:     SNESSetFunction(snes,r,FormFunction1_block,NULL);

428:     /*
429:     Set Jacobian matrix data structure and Jacobian evaluation routine
430:     */
431:     SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
432:   } else {
433:     SNESSetFunction(snes,r,FormFunction2_block,NULL);
434:     SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
435:   }

437:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
438:   Customize nonlinear solver; set runtime options
439:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

451:   /*
452:   Set SNES/KSP/KSP/PC runtime options, e.g.,
453:   -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
454:   These options will override those specified above as long as
455:   SNESSetFromOptions() is called _after_ any other customization
456:   routines.
457:   */
458:   SNESSetFromOptions(snes);

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

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

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

497:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
498:   Free work space.  All PETSc objects should be destroyed when they
499:   are no longer needed.
500:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
501:   VecDestroy(&x)); PetscCall(VecDestroy(&r);
502:   MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
503:   return 0;
504: }

506: /* ------------------------------------------------------------------- */
507: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
508: {
509:   Vec            *xx, *ff, x1,x2, f1,f2;
510:   PetscScalar    ff_0, ff_1;
511:   PetscScalar    xx_0, xx_1;
512:   PetscInt       index,nb;

515:   /* get blocks for function */
516:   VecNestGetSubVecs(f, &nb, &ff);
517:   f1   = ff[0];  f2 = ff[1];

519:   /* get blocks for solution */
520:   VecNestGetSubVecs(x, &nb, &xx);
521:   x1   = xx[0];  x2 = xx[1];

523:   /* get solution values */
524:   index = 0;
525:   VecGetValues(x1,1, &index, &xx_0);
526:   VecGetValues(x2,1, &index, &xx_1);

528:   /* Compute function */
529:   ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
530:   ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;

532:   /* set function values */
533:   VecSetValue(f1, index, ff_0, INSERT_VALUES);

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

537:   VecAssemblyBegin(f);
538:   VecAssemblyEnd(f);
539:   return 0;
540: }

542: /* ------------------------------------------------------------------- */
543: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
544: {
545:   Vec            *xx, x1,x2;
546:   PetscScalar    xx_0, xx_1;
547:   PetscInt       index,nb;
548:   PetscScalar    A_00, A_01, A_10, A_11;
549:   Mat            j11, j12, j21, j22;
550:   Mat            **mats;

553:   /* get blocks for solution */
554:   VecNestGetSubVecs(x, &nb, &xx);
555:   x1   = xx[0];  x2 = xx[1];

557:   /* get solution values */
558:   index = 0;
559:   VecGetValues(x1,1, &index, &xx_0);
560:   VecGetValues(x2,1, &index, &xx_1);

562:   /* get block matrices */
563:   MatNestGetSubMats(jac,NULL,NULL,&mats);
564:   j11  = mats[0][0];
565:   j12  = mats[0][1];
566:   j21  = mats[1][0];
567:   j22  = mats[1][1];

569:   /* compute jacobian entries */
570:   A_00 = 2.0*xx_0 + xx_1;
571:   A_01 = xx_0;
572:   A_10 = xx_1;
573:   A_11 = xx_0 + 2.0*xx_1;

575:   /* set jacobian values */
576:   MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
577:   MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
578:   MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
579:   MatSetValue(j22, 0,0, A_11, INSERT_VALUES);

581:   /* Assemble sub matrix */
582:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
583:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
584:   return 0;
585: }

587: /* ------------------------------------------------------------------- */
588: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
589: {
590:   PetscScalar       *ff;
591:   const PetscScalar *xx;

594:   /*
595:   Get pointers to vector data.
596:   - For default PETSc vectors, VecGetArray() returns a pointer to
597:   the data array.  Otherwise, the routine is implementation dependent.
598:   - You MUST call VecRestoreArray() when you no longer need access to
599:   the array.
600:   */
601:   VecGetArrayRead(x,&xx);
602:   VecGetArray(f,&ff);

604:   /*
605:   Compute function
606:   */
607:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
608:   ff[1] = xx[1];

610:   /*
611:   Restore vectors
612:   */
613:   VecRestoreArrayRead(x,&xx);
614:   VecRestoreArray(f,&ff);
615:   return 0;
616: }

618: /* ------------------------------------------------------------------- */
619: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
620: {
621:   const PetscScalar *xx;
622:   PetscScalar       A[4];
623:   PetscInt          idx[2] = {0,1};

626:   /*
627:   Get pointer to vector data
628:   */
629:   VecGetArrayRead(x,&xx);

631:   /*
632:   Compute Jacobian entries and insert into matrix.
633:   - Since this is such a small problem, we set all entries for
634:   the matrix at once.
635:   */
636:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
637:   A[2]  = 0.0;                     A[3] = 1.0;
638:   MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);

640:   /*
641:   Restore vector
642:   */
643:   VecRestoreArrayRead(x,&xx);

645:   /*
646:   Assemble matrix
647:   */
648:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
649:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
650:   return 0;
651: }

653: int main(int argc,char **argv)
654: {
655:   PetscMPIInt    size;

657:   PetscInitialize(&argc,&argv,(char*)0,help);
658:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
660:   assembled_system();
661:   block_system();
662:   PetscFinalize();
663:   return 0;
664: }

666: /*TEST

668:    test:
669:       args: -snes_monitor_short
670:       requires: !single

672: TEST*/