Actual source code: ex17.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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,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*);


 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(NULL,NULL,"-hard",&flg);
100:   if (!flg) {
101:     /*
102:     Set function evaluation routine and vector.
103:     */
104:     SNESSetFunction(snes,r,FormFunction1,NULL);

106:     /*
107:     Set Jacobian matrix data structure and Jacobian evaluation routine
108:     */
109:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
110:   } else {
111:     SNESSetFunction(snes,r,FormFunction2,NULL);
112:     SNESSetJacobian(snes,J,J,FormJacobian2,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,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: }

176: /* ------------------------------------------------------------------- */
179: /*
180: FormFunction1 - Evaluates nonlinear function, F(x).

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

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

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

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


214:   /*
215:   Restore vectors
216:   */
217:   VecRestoreArrayRead(x,&xx);
218:   VecRestoreArray(f,&ff);
219:   return(0);
220: }
221: /* ------------------------------------------------------------------- */
224: /*
225: FormJacobian1 - Evaluates Jacobian matrix.

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

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

245:   /*
246:   Get pointer to vector data
247:   */
248:   VecGetArrayRead(x,&xx);

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

259:   /*
260:   Restore vector
261:   */
262:   VecRestoreArrayRead(x,&xx);

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


273: /* ------------------------------------------------------------------- */
276: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
277: {
278:   PetscErrorCode    ierr;
279:   const PetscScalar *xx;
280:   PetscScalar       *ff;

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

293:   /*
294:   Compute function
295:   */
296:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
297:   ff[1] = xx[1];

299:   /*
300:   Restore vectors
301:   */
302:   VecRestoreArrayRead(x,&xx);
303:   VecRestoreArray(f,&ff);
304:   return(0);
305: }

307: /* ------------------------------------------------------------------- */
310: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
311: {
312:   const PetscScalar *xx;
313:   PetscScalar       A[4];
314:   PetscErrorCode    ierr;
315:   PetscInt          idx[2] = {0,1};

318:   /*
319:   Get pointer to vector data
320:   */
321:   VecGetArrayRead(x,&xx);

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

332:   /*
333:   Restore vector
334:   */
335:   VecRestoreArrayRead(x,&xx);

337:   /*
338:   Assemble matrix
339:   */
340:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
341:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
342:   return(0);
343: }

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

359:   Mat j11, j12, j21, j22;
360:   Vec x1, x2, r1, r2;
361:   Vec bv;
362:   Vec bx[2];
363:   Mat bA[2][2];

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

368:   SNESCreate(PETSC_COMM_WORLD,&snes);

370:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371:   Create matrix and vector data structures; set corresponding routines
372:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

374:   /*
375:   Create sub vectors for solution and nonlinear function
376:   */
377:   VecCreateSeq(PETSC_COMM_SELF,1,&x1);
378:   VecDuplicate(x1,&r1);

380:   VecCreateSeq(PETSC_COMM_SELF,1,&x2);
381:   VecDuplicate(x2,&r2);

383:   /*
384:   Create the block vectors
385:   */
386:   bx[0] = x1;
387:   bx[1] = x2;
388:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
389:   VecAssemblyBegin(x);
390:   VecAssemblyEnd(x);
391:   VecDestroy(&x1);
392:   VecDestroy(&x2);

394:   bx[0] = r1;
395:   bx[1] = r2;
396:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
397:   VecDestroy(&r1);
398:   VecDestroy(&r2);
399:   VecAssemblyBegin(r);
400:   VecAssemblyEnd(r);

402:   /*
403:   Create sub Jacobian matrix data structure
404:   */
405:   MatCreate(PETSC_COMM_WORLD, &j11);
406:   MatSetSizes(j11, 1, 1, 1, 1);
407:   MatSetType(j11, MATSEQAIJ);
408:   MatSetUp(j11);

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

415:   MatCreate(PETSC_COMM_WORLD, &j21);
416:   MatSetSizes(j21, 1, 1, 1, 1);
417:   MatSetType(j21, MATSEQAIJ);
418:   MatSetUp(j21);

420:   MatCreate(PETSC_COMM_WORLD, &j22);
421:   MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
422:   MatSetType(j22, MATSEQAIJ);
423:   MatSetUp(j22);
424:   /*
425:   Create block Jacobian matrix data structure
426:   */
427:   bA[0][0] = j11;
428:   bA[0][1] = j12;
429:   bA[1][0] = j21;
430:   bA[1][1] = j22;

432:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
433:   MatSetUp(J);
434:   MatNestSetVecType(J,VECNEST);
435:   MatDestroy(&j11);
436:   MatDestroy(&j12);
437:   MatDestroy(&j21);
438:   MatDestroy(&j22);

440:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
441:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

443:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
444:   if (!flg) {
445:     /*
446:     Set function evaluation routine and vector.
447:     */
448:     SNESSetFunction(snes,r,FormFunction1_block,NULL);

450:     /*
451:     Set Jacobian matrix data structure and Jacobian evaluation routine
452:     */
453:     SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
454:   } else {
455:     SNESSetFunction(snes,r,FormFunction2_block,NULL);
456:     SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
457:   }

459:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460:   Customize nonlinear solver; set runtime options
461:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

519:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
520:   Free work space.  All PETSc objects should be destroyed when they
521:   are no longer needed.
522:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
523:   VecDestroy(&x); VecDestroy(&r);
524:   MatDestroy(&J); SNESDestroy(&snes);
525:   return(0);
526: }

528: /* ------------------------------------------------------------------- */
531: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
532: {
533:   Vec            *xx, *ff, x1,x2, f1,f2;
534:   PetscScalar    ff_0, ff_1;
535:   PetscScalar    xx_0, xx_1;
536:   PetscInt       index,nb;

540:   /* get blocks for function */
541:   VecNestGetSubVecs(f, &nb, &ff);
542:   f1   = ff[0];  f2 = ff[1];

544:   /* get blocks for solution */
545:   VecNestGetSubVecs(x, &nb, &xx);
546:   x1   = xx[0];  x2 = xx[1];

548:   /* get solution values */
549:   index = 0;
550:   VecGetValues(x1,1, &index, &xx_0);
551:   VecGetValues(x2,1, &index, &xx_1);

553:   /* Compute function */
554:   ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
555:   ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;

557:   /* set function values */
558:   VecSetValue(f1, index, ff_0, INSERT_VALUES);

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

562:   VecAssemblyBegin(f);
563:   VecAssemblyEnd(f);
564:   return(0);
565: }

567: /* ------------------------------------------------------------------- */
570: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
571: {
572:   Vec            *xx, x1,x2;
573:   PetscScalar    xx_0, xx_1;
574:   PetscInt       index,nb;
575:   PetscScalar    A_00, A_01, A_10, A_11;
576:   Mat            j11, j12, j21, j22;
577:   Mat            **mats;

581:   /* get blocks for solution */
582:   VecNestGetSubVecs(x, &nb, &xx);
583:   x1   = xx[0];  x2 = xx[1];

585:   /* get solution values */
586:   index = 0;
587:   VecGetValues(x1,1, &index, &xx_0);
588:   VecGetValues(x2,1, &index, &xx_1);

590:   /* get block matrices */
591:   MatNestGetSubMats(jac,NULL,NULL,&mats);
592:   j11  = mats[0][0];
593:   j12  = mats[0][1];
594:   j21  = mats[1][0];
595:   j22  = mats[1][1];

597:   /* compute jacobian entries */
598:   A_00 = 2.0*xx_0 + xx_1;
599:   A_01 = xx_0;
600:   A_10 = xx_1;
601:   A_11 = xx_0 + 2.0*xx_1;

603:   /* set jacobian values */
604:   MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
605:   MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
606:   MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
607:   MatSetValue(j22, 0,0, A_11, INSERT_VALUES);

609:   /* Assemble sub matrix */
610:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
611:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
612:   return(0);
613: }

615: /* ------------------------------------------------------------------- */
618: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
619: {
620:   PetscErrorCode    ierr;
621:   PetscScalar       *ff;
622:   const PetscScalar *xx;

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

635:   /*
636:   Compute function
637:   */
638:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
639:   ff[1] = xx[1];

641:   /*
642:   Restore vectors
643:   */
644:   VecRestoreArrayRead(x,&xx);
645:   VecRestoreArray(f,&ff);
646:   return(0);
647: }

649: /* ------------------------------------------------------------------- */
652: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
653: {
654:   const PetscScalar *xx;
655:   PetscScalar       A[4];
656:   PetscErrorCode    ierr;
657:   PetscInt          idx[2] = {0,1};

660:   /*
661:   Get pointer to vector data
662:   */
663:   VecGetArrayRead(x,&xx);

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

674:   /*
675:   Restore vector
676:   */
677:   VecRestoreArrayRead(x,&xx);

679:   /*
680:   Assemble matrix
681:   */
682:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
683:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
684:   return(0);
685: }

689: int main(int argc,char **argv)
690: {
691:   PetscMPIInt    size;

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

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

699:   assembled_system();

701:   block_system();

703:   PetscFinalize();
704:   return 0;
705: }