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: /*
  5: Include "petscsnes.h" so that we can use SNES solvers.  Note that this
  6: file automatically includes:
  7: petscsys.h       - base PETSc routines   petscvec.h - vectors
  8: petscsys.h    - system routines       petscmat.h - matrices
  9: petscis.h     - index sets            petscksp.h - Krylov subspace methods
 10: petscviewer.h - viewers               petscpc.h  - preconditioners
 11: petscksp.h   - linear solvers
 12: */
 13: #include <petscsnes.h>

 15: /*
 16: This example is block version of the test found at
 17:   ${PETSC_DIR}/src/snes/tutorials/ex1.c
 18: In this test we replace the Jacobian systems
 19:   [J]{x} = {F}
 20: where

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

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

 29: [J] = ([j00], [j01]),  {x} = ({x0}, {x1})^T, {F} = ({f0}, {f1})^T
 30:       ([j10], [j11])
 31: where
 32: [j00], [j01], [j10], [j11] \in \mathbb^{1 \times 1}
 33: {x0}, {x1}, {f0}, {f1} \in \mathbb^{1 \times 1}

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

 40: /*
 41: User-defined routines
 42: */
 43: static PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
 44: static PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 45: static PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
 46: static PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
 47: static PetscErrorCode FormJacobian1_block(SNES,Vec,Mat,Mat,void*);
 48: static PetscErrorCode FormFunction1_block(SNES,Vec,Vec,void*);
 49: static PetscErrorCode FormJacobian2_block(SNES,Vec,Mat,Mat,void*);
 50: static PetscErrorCode FormFunction2_block(SNES,Vec,Vec,void*);

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

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

 66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:   Create nonlinear solver context
 68:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 70:   SNESCreate(PETSC_COMM_WORLD,&snes);

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 73:   Create matrix and vector data structures; set corresponding routines
 74:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 76:   /*
 77:   Create vectors for solution and nonlinear function
 78:   */
 79:   VecCreateSeq(PETSC_COMM_SELF,2,&x);
 80:   VecDuplicate(x,&r);

 82:   /*
 83:   Create Jacobian matrix data structure
 84:   */
 85:   MatCreate(PETSC_COMM_SELF,&J);
 86:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 87:   MatSetFromOptions(J);
 88:   MatSetUp(J);

 90:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
 91:   if (!flg) {
 92:     /*
 93:     Set function evaluation routine and vector.
 94:     */
 95:     SNESSetFunction(snes,r,FormFunction1,NULL);

 97:     /*
 98:     Set Jacobian matrix data structure and Jacobian evaluation routine
 99:     */
100:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
101:   } else {
102:     SNESSetFunction(snes,r,FormFunction2,NULL);
103:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
104:   }

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:   Customize nonlinear solver; set runtime options
108:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

146:   SNESSolve(snes,NULL,x);
147:   SNESGetIterationNumber(snes,&its);
148:   if (flg) {
149:     Vec f;
150:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
151:     SNESGetFunction(snes,&f,0,0);
152:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
153:   }
154:   PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);

156:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157:   Free work space.  All PETSc objects should be destroyed when they
158:   are no longer needed.
159:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

161:   VecDestroy(&x)); PetscCall(VecDestroy(&r);
162:   MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
163:   return 0;
164: }

166: /* ------------------------------------------------------------------- */
167: /*
168: FormFunction1 - Evaluates nonlinear function, F(x).

170: Input Parameters:
171: .  snes - the SNES context
172: .  x - input vector
173: .  dummy - optional user-defined context (not used here)

175: Output Parameter:
176: .  f - function vector
177: */
178: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
179: {
180:   const PetscScalar *xx;
181:   PetscScalar       *ff;

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

194:   /*
195:   Compute function
196:   */
197:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
198:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

200:   /*
201:   Restore vectors
202:   */
203:   VecRestoreArrayRead(x,&xx);
204:   VecRestoreArray(f,&ff);
205:   return 0;
206: }
207: /* ------------------------------------------------------------------- */
208: /*
209: FormJacobian1 - Evaluates Jacobian matrix.

211: Input Parameters:
212: .  snes - the SNES context
213: .  x - input vector
214: .  dummy - optional user-defined context (not used here)

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

228:   /*
229:   Get pointer to vector data
230:   */
231:   VecGetArrayRead(x,&xx);

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

242:   /*
243:   Restore vector
244:   */
245:   VecRestoreArrayRead(x,&xx);

247:   /*
248:   Assemble matrix
249:   */
250:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
251:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
252:   return 0;
253: }

255: /* ------------------------------------------------------------------- */
256: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
257: {
258:   const PetscScalar *xx;
259:   PetscScalar       *ff;

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

272:   /*
273:   Compute function
274:   */
275:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
276:   ff[1] = xx[1];

278:   /*
279:   Restore vectors
280:   */
281:   VecRestoreArrayRead(x,&xx);
282:   VecRestoreArray(f,&ff);
283:   return 0;
284: }

286: /* ------------------------------------------------------------------- */
287: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
288: {
289:   const PetscScalar *xx;
290:   PetscScalar       A[4];
291:   PetscInt          idx[2] = {0,1};

294:   /*
295:   Get pointer to vector data
296:   */
297:   VecGetArrayRead(x,&xx);

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

308:   /*
309:   Restore vector
310:   */
311:   VecRestoreArrayRead(x,&xx);

313:   /*
314:   Assemble matrix
315:   */
316:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
317:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
318:   return 0;
319: }

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

332:   Mat            j11, j12, j21, j22;
333:   Vec            x1, x2, r1, r2;
334:   Vec            bv;
335:   Vec            bx[2];
336:   Mat            bA[2][2];

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

341:   SNESCreate(PETSC_COMM_WORLD,&snes);

343:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
344:   Create matrix and vector data structures; set corresponding routines
345:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

347:   /*
348:   Create sub vectors for solution and nonlinear function
349:   */
350:   VecCreateSeq(PETSC_COMM_SELF,1,&x1);
351:   VecDuplicate(x1,&r1);

353:   VecCreateSeq(PETSC_COMM_SELF,1,&x2);
354:   VecDuplicate(x2,&r2);

356:   /*
357:   Create the block vectors
358:   */
359:   bx[0] = x1;
360:   bx[1] = x2;
361:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
362:   VecAssemblyBegin(x);
363:   VecAssemblyEnd(x);
364:   VecDestroy(&x1);
365:   VecDestroy(&x2);

367:   bx[0] = r1;
368:   bx[1] = r2;
369:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
370:   VecDestroy(&r1);
371:   VecDestroy(&r2);
372:   VecAssemblyBegin(r);
373:   VecAssemblyEnd(r);

375:   /*
376:   Create sub Jacobian matrix data structure
377:   */
378:   MatCreate(PETSC_COMM_WORLD, &j11);
379:   MatSetSizes(j11, 1, 1, 1, 1);
380:   MatSetType(j11, MATSEQAIJ);
381:   MatSetUp(j11);

383:   MatCreate(PETSC_COMM_WORLD, &j12);
384:   MatSetSizes(j12, 1, 1, 1, 1);
385:   MatSetType(j12, MATSEQAIJ);
386:   MatSetUp(j12);

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

393:   MatCreate(PETSC_COMM_WORLD, &j22);
394:   MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
395:   MatSetType(j22, MATSEQAIJ);
396:   MatSetUp(j22);
397:   /*
398:   Create block Jacobian matrix data structure
399:   */
400:   bA[0][0] = j11;
401:   bA[0][1] = j12;
402:   bA[1][0] = j21;
403:   bA[1][1] = j22;

405:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
406:   MatSetUp(J);
407:   MatNestSetVecType(J,VECNEST);
408:   MatDestroy(&j11);
409:   MatDestroy(&j12);
410:   MatDestroy(&j21);
411:   MatDestroy(&j22);

413:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
414:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

416:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
417:   if (!flg) {
418:     /*
419:     Set function evaluation routine and vector.
420:     */
421:     SNESSetFunction(snes,r,FormFunction1_block,NULL);

423:     /*
424:     Set Jacobian matrix data structure and Jacobian evaluation routine
425:     */
426:     SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
427:   } else {
428:     SNESSetFunction(snes,r,FormFunction2_block,NULL);
429:     SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
430:   }

432:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
433:   Customize nonlinear solver; set runtime options
434:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

492:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
493:   Free work space.  All PETSc objects should be destroyed when they
494:   are no longer needed.
495:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
496:   VecDestroy(&x)); PetscCall(VecDestroy(&r);
497:   MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
498:   return 0;
499: }

501: /* ------------------------------------------------------------------- */
502: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
503: {
504:   Vec            *xx, *ff, x1,x2, f1,f2;
505:   PetscScalar    ff_0, ff_1;
506:   PetscScalar    xx_0, xx_1;
507:   PetscInt       index,nb;

510:   /* get blocks for function */
511:   VecNestGetSubVecs(f, &nb, &ff);
512:   f1   = ff[0];  f2 = ff[1];

514:   /* get blocks for solution */
515:   VecNestGetSubVecs(x, &nb, &xx);
516:   x1   = xx[0];  x2 = xx[1];

518:   /* get solution values */
519:   index = 0;
520:   VecGetValues(x1,1, &index, &xx_0);
521:   VecGetValues(x2,1, &index, &xx_1);

523:   /* Compute function */
524:   ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
525:   ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;

527:   /* set function values */
528:   VecSetValue(f1, index, ff_0, INSERT_VALUES);

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

532:   VecAssemblyBegin(f);
533:   VecAssemblyEnd(f);
534:   return 0;
535: }

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

548:   /* get blocks for solution */
549:   VecNestGetSubVecs(x, &nb, &xx);
550:   x1   = xx[0];  x2 = xx[1];

552:   /* get solution values */
553:   index = 0;
554:   VecGetValues(x1,1, &index, &xx_0);
555:   VecGetValues(x2,1, &index, &xx_1);

557:   /* get block matrices */
558:   MatNestGetSubMats(jac,NULL,NULL,&mats);
559:   j11  = mats[0][0];
560:   j12  = mats[0][1];
561:   j21  = mats[1][0];
562:   j22  = mats[1][1];

564:   /* compute jacobian entries */
565:   A_00 = 2.0*xx_0 + xx_1;
566:   A_01 = xx_0;
567:   A_10 = xx_1;
568:   A_11 = xx_0 + 2.0*xx_1;

570:   /* set jacobian values */
571:   MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
572:   MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
573:   MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
574:   MatSetValue(j22, 0,0, A_11, INSERT_VALUES);

576:   /* Assemble sub matrix */
577:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
578:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
579:   return 0;
580: }

582: /* ------------------------------------------------------------------- */
583: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
584: {
585:   PetscScalar       *ff;
586:   const PetscScalar *xx;

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

599:   /*
600:   Compute function
601:   */
602:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
603:   ff[1] = xx[1];

605:   /*
606:   Restore vectors
607:   */
608:   VecRestoreArrayRead(x,&xx);
609:   VecRestoreArray(f,&ff);
610:   return 0;
611: }

613: /* ------------------------------------------------------------------- */
614: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
615: {
616:   const PetscScalar *xx;
617:   PetscScalar       A[4];
618:   PetscInt          idx[2] = {0,1};

621:   /*
622:   Get pointer to vector data
623:   */
624:   VecGetArrayRead(x,&xx);

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

635:   /*
636:   Restore vector
637:   */
638:   VecRestoreArrayRead(x,&xx);

640:   /*
641:   Assemble matrix
642:   */
643:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
644:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
645:   return 0;
646: }

648: int main(int argc,char **argv)
649: {
650:   PetscMPIInt    size;

652:   PetscInitialize(&argc,&argv,(char*)0,help);
653:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
655:   assembled_system();
656:   block_system();
657:   PetscFinalize();
658:   return 0;
659: }

661: /*TEST

663:    test:
664:       args: -snes_monitor_short
665:       requires: !single

667: TEST*/