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;

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

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

 70:   PetscCall(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:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, 2, &x));
 80:   PetscCall(VecDuplicate(x, &r));

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

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

 97:     /*
 98:     Set Jacobian matrix data structure and Jacobian evaluation routine
 99:     */
100:     PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
101:   } else {
102:     PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL));
103:     PetscCall(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:   PetscCall(SNESGetKSP(snes, &ksp));
116:   PetscCall(KSPGetPC(ksp, &pc));
117:   PetscCall(PCSetType(pc, PCNONE));
118:   PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_CURRENT, PETSC_CURRENT, 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:   PetscCall(SNESSetFromOptions(snes));

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

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

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

162:   PetscCall(VecDestroy(&x));
163:   PetscCall(VecDestroy(&r));
164:   PetscCall(MatDestroy(&J));
165:   PetscCall(SNESDestroy(&snes));
166:   PetscFunctionReturn(PETSC_SUCCESS);
167: }

169: /*
170: FormFunction1 - Evaluates nonlinear function, F(x).

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

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

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

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

202:   /*
203:   Restore vectors
204:   */
205:   PetscCall(VecRestoreArrayRead(x, &xx));
206:   PetscCall(VecRestoreArray(f, &ff));
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: /*
211: FormJacobian1 - Evaluates Jacobian matrix.

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

218: Output Parameters:
219: .  jac - Jacobian matrix
220: .  B - optionally different preconditioning matrix

222: */
223: static PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
224: {
225:   const PetscScalar *xx;
226:   PetscScalar        A[4];
227:   PetscInt           idx[2] = {0, 1};

229:   PetscFunctionBeginUser;
230:   /*
231:   Get pointer to vector data
232:   */
233:   PetscCall(VecGetArrayRead(x, &xx));

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

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

251:   /*
252:   Assemble matrix
253:   */
254:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
255:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

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

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

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

281:   /*
282:   Restore vectors
283:   */
284:   PetscCall(VecRestoreArrayRead(x, &xx));
285:   PetscCall(VecRestoreArray(f, &ff));
286:   PetscFunctionReturn(PETSC_SUCCESS);
287: }

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

295:   PetscFunctionBeginUser;
296:   /*
297:   Get pointer to vector data
298:   */
299:   PetscCall(VecGetArrayRead(x, &xx));

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

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

317:   /*
318:   Assemble matrix
319:   */
320:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
321:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

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

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

342:   PetscFunctionBeginUser;
343:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n"));

345:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

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

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

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

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

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

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

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

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

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

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

417:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
418:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));

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

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

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

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

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

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

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

494:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));

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

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;

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

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

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

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

534:   /* set function values */
535:   PetscCall(VecSetValue(f1, index, ff_0, INSERT_VALUES));

537:   PetscCall(VecSetValue(f2, index, ff_1, INSERT_VALUES));

539:   PetscCall(VecAssemblyBegin(f));
540:   PetscCall(VecAssemblyEnd(f));
541:   PetscFunctionReturn(PETSC_SUCCESS);
542: }

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

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

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

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

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

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

583:   /* Assemble sub matrix */
584:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
585:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
586:   PetscFunctionReturn(PETSC_SUCCESS);
587: }

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

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

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

611:   /*
612:   Restore vectors
613:   */
614:   PetscCall(VecRestoreArrayRead(x, &xx));
615:   PetscCall(VecRestoreArray(f, &ff));
616:   PetscFunctionReturn(PETSC_SUCCESS);
617: }

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};

625:   PetscFunctionBeginUser;
626:   /*
627:   Get pointer to vector data
628:   */
629:   PetscCall(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;
637:   A[1] = 0.0;
638:   A[2] = 0.0;
639:   A[3] = 1.0;
640:   PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));

642:   /*
643:   Restore vector
644:   */
645:   PetscCall(VecRestoreArrayRead(x, &xx));

647:   /*
648:   Assemble matrix
649:   */
650:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
651:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
652:   PetscFunctionReturn(PETSC_SUCCESS);
653: }

655: int main(int argc, char **argv)
656: {
657:   PetscMPIInt size;

659:   PetscFunctionBeginUser;
660:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
661:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
662:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
663:   PetscCall(assembled_system());
664:   PetscCall(block_system());
665:   PetscCall(PetscFinalize());
666:   return 0;
667: }

669: /*TEST

671:    test:
672:       args: -snes_monitor_short
673:       requires: !single

675: TEST*/