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) PetscCall(VecSet(x, pfive));
133: else {
134: PetscCall(VecGetArray(x, &xx));
135: xx[0] = 2.0;
136: xx[1] = 3.0;
137: PetscCall(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: PetscCall(SNESSolve(snes, NULL, x));
147: PetscCall(SNESGetIterationNumber(snes, &its));
148: if (flg) {
149: Vec f;
150: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
151: PetscCall(SNESGetFunction(snes, &f, 0, 0));
152: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
153: }
154: PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Free work space. All PETSc objects should be destroyed when they
158: are no longer needed.
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: PetscCall(VecDestroy(&x));
162: PetscCall(VecDestroy(&r));
163: PetscCall(MatDestroy(&J));
164: PetscCall(SNESDestroy(&snes));
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*
169: FormFunction1 - Evaluates nonlinear function, F(x).
171: Input Parameters:
172: . snes - the SNES context
173: . x - input vector
174: . dummy - optional user-defined context (not used here)
176: Output Parameter:
177: . f - function vector
178: */
179: static PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *dummy)
180: {
181: const PetscScalar *xx;
182: PetscScalar *ff;
184: PetscFunctionBeginUser;
185: /*
186: Get pointers to vector data.
187: - For default PETSc vectors, VecGetArray() returns a pointer to
188: the data array. Otherwise, the routine is implementation dependent.
189: - You MUST call VecRestoreArray() when you no longer need access to
190: the array.
191: */
192: PetscCall(VecGetArrayRead(x, &xx));
193: PetscCall(VecGetArray(f, &ff));
195: /*
196: Compute function
197: */
198: ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0;
199: ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0;
201: /*
202: Restore vectors
203: */
204: PetscCall(VecRestoreArrayRead(x, &xx));
205: PetscCall(VecRestoreArray(f, &ff));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*
210: FormJacobian1 - Evaluates Jacobian matrix.
212: Input Parameters:
213: . snes - the SNES context
214: . x - input vector
215: . dummy - optional user-defined context (not used here)
217: Output Parameters:
218: . jac - Jacobian matrix
219: . B - optionally different matrix used to construct the preconditioner
221: */
222: static PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
223: {
224: const PetscScalar *xx;
225: PetscScalar A[4];
226: PetscInt idx[2] = {0, 1};
228: PetscFunctionBeginUser;
229: /*
230: Get pointer to vector data
231: */
232: PetscCall(VecGetArrayRead(x, &xx));
234: /*
235: Compute Jacobian entries and insert into matrix.
236: - Since this is such a small problem, we set all entries for
237: the matrix at once.
238: */
239: A[0] = 2.0 * xx[0] + xx[1];
240: A[1] = xx[0];
241: A[2] = xx[1];
242: A[3] = xx[0] + 2.0 * xx[1];
243: PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
245: /*
246: Restore vector
247: */
248: PetscCall(VecRestoreArrayRead(x, &xx));
250: /*
251: Assemble matrix
252: */
253: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
254: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: static PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
259: {
260: const PetscScalar *xx;
261: PetscScalar *ff;
263: PetscFunctionBeginUser;
264: /*
265: Get pointers to vector data.
266: - For default PETSc vectors, VecGetArray() returns a pointer to
267: the data array. Otherwise, the routine is implementation dependent.
268: - You MUST call VecRestoreArray() when you no longer need access to
269: the array.
270: */
271: PetscCall(VecGetArrayRead(x, &xx));
272: PetscCall(VecGetArray(f, &ff));
274: /*
275: Compute function
276: */
277: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
278: ff[1] = xx[1];
280: /*
281: Restore vectors
282: */
283: PetscCall(VecRestoreArrayRead(x, &xx));
284: PetscCall(VecRestoreArray(f, &ff));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: static PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
289: {
290: const PetscScalar *xx;
291: PetscScalar A[4];
292: PetscInt idx[2] = {0, 1};
294: PetscFunctionBeginUser;
295: /*
296: Get pointer to vector data
297: */
298: PetscCall(VecGetArrayRead(x, &xx));
300: /*
301: Compute Jacobian entries and insert into matrix.
302: - Since this is such a small problem, we set all entries for
303: the matrix at once.
304: */
305: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
306: A[1] = 0.0;
307: A[2] = 0.0;
308: A[3] = 1.0;
309: PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
311: /*
312: Restore vector
313: */
314: PetscCall(VecRestoreArrayRead(x, &xx));
316: /*
317: Assemble matrix
318: */
319: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
320: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: static PetscErrorCode block_system(void)
325: {
326: SNES snes; /* nonlinear solver context */
327: KSP ksp; /* linear solver context */
328: PC pc; /* preconditioner context */
329: Vec x, r; /* solution, residual vectors */
330: Mat J; /* Jacobian matrix */
331: PetscInt its;
332: PetscScalar pfive = .5;
333: PetscBool flg;
335: Mat j11, j12, j21, j22;
336: Vec x1, x2, r1, r2;
337: Vec bv;
338: Vec bx[2];
339: Mat bA[2][2];
341: PetscFunctionBeginUser;
342: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n"));
344: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
346: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
347: Create matrix and vector data structures; set corresponding routines
348: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350: /*
351: Create sub vectors for solution and nonlinear function
352: */
353: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &x1));
354: PetscCall(VecDuplicate(x1, &r1));
356: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &x2));
357: PetscCall(VecDuplicate(x2, &r2));
359: /*
360: Create the block vectors
361: */
362: bx[0] = x1;
363: bx[1] = x2;
364: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &x));
365: PetscCall(VecAssemblyBegin(x));
366: PetscCall(VecAssemblyEnd(x));
367: PetscCall(VecDestroy(&x1));
368: PetscCall(VecDestroy(&x2));
370: bx[0] = r1;
371: bx[1] = r2;
372: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &r));
373: PetscCall(VecDestroy(&r1));
374: PetscCall(VecDestroy(&r2));
375: PetscCall(VecAssemblyBegin(r));
376: PetscCall(VecAssemblyEnd(r));
378: /*
379: Create sub Jacobian matrix data structure
380: */
381: PetscCall(MatCreate(PETSC_COMM_WORLD, &j11));
382: PetscCall(MatSetSizes(j11, 1, 1, 1, 1));
383: PetscCall(MatSetType(j11, MATSEQAIJ));
384: PetscCall(MatSetUp(j11));
386: PetscCall(MatCreate(PETSC_COMM_WORLD, &j12));
387: PetscCall(MatSetSizes(j12, 1, 1, 1, 1));
388: PetscCall(MatSetType(j12, MATSEQAIJ));
389: PetscCall(MatSetUp(j12));
391: PetscCall(MatCreate(PETSC_COMM_WORLD, &j21));
392: PetscCall(MatSetSizes(j21, 1, 1, 1, 1));
393: PetscCall(MatSetType(j21, MATSEQAIJ));
394: PetscCall(MatSetUp(j21));
396: PetscCall(MatCreate(PETSC_COMM_WORLD, &j22));
397: PetscCall(MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
398: PetscCall(MatSetType(j22, MATSEQAIJ));
399: PetscCall(MatSetUp(j22));
400: /*
401: Create block Jacobian matrix data structure
402: */
403: bA[0][0] = j11;
404: bA[0][1] = j12;
405: bA[1][0] = j21;
406: bA[1][1] = j22;
408: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &bA[0][0], &J));
409: PetscCall(MatSetUp(J));
410: PetscCall(MatNestSetVecType(J, VECNEST));
411: PetscCall(MatDestroy(&j11));
412: PetscCall(MatDestroy(&j12));
413: PetscCall(MatDestroy(&j21));
414: PetscCall(MatDestroy(&j22));
416: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
417: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
419: PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg));
420: if (!flg) {
421: /*
422: Set function evaluation routine and vector.
423: */
424: PetscCall(SNESSetFunction(snes, r, FormFunction1_block, NULL));
426: /*
427: Set Jacobian matrix data structure and Jacobian evaluation routine
428: */
429: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1_block, NULL));
430: } else {
431: PetscCall(SNESSetFunction(snes, r, FormFunction2_block, NULL));
432: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2_block, NULL));
433: }
435: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436: Customize nonlinear solver; set runtime options
437: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
439: /*
440: Set linear solver defaults for this problem. By extracting the
441: KSP, KSP, and PC contexts from the SNES context, we can then
442: directly call any KSP, KSP, and PC routines to set various options.
443: */
444: PetscCall(SNESGetKSP(snes, &ksp));
445: PetscCall(KSPGetPC(ksp, &pc));
446: PetscCall(PCSetType(pc, PCNONE));
447: PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_CURRENT, PETSC_CURRENT, 20));
449: /*
450: Set SNES/KSP/KSP/PC runtime options, e.g.,
451: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
452: These options will override those specified above as long as
453: SNESSetFromOptions() is called _after_ any other customization
454: routines.
455: */
456: PetscCall(SNESSetFromOptions(snes));
458: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
459: Evaluate initial guess; then solve nonlinear system
460: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
461: if (!flg) PetscCall(VecSet(x, pfive));
462: else {
463: Vec *vecs;
464: PetscCall(VecNestGetSubVecs(x, NULL, &vecs));
465: bv = vecs[0];
466: /* PetscCall(VecBlockGetSubVec(x, 0, &bv)); */
467: PetscCall(VecSetValue(bv, 0, 2.0, INSERT_VALUES)); /* xx[0] = 2.0; */
468: PetscCall(VecAssemblyBegin(bv));
469: PetscCall(VecAssemblyEnd(bv));
471: /* PetscCall(VecBlockGetSubVec(x, 1, &bv)); */
472: bv = vecs[1];
473: PetscCall(VecSetValue(bv, 0, 3.0, INSERT_VALUES)); /* xx[1] = 3.0; */
474: PetscCall(VecAssemblyBegin(bv));
475: PetscCall(VecAssemblyEnd(bv));
476: }
477: /*
478: Note: The user should initialize the vector, x, with the initial guess
479: for the nonlinear solver prior to calling SNESSolve(). In particular,
480: to employ an initial guess of zero, the user should explicitly set
481: this vector to zero by calling VecSet().
482: */
483: PetscCall(SNESSolve(snes, NULL, x));
484: PetscCall(SNESGetIterationNumber(snes, &its));
485: if (flg) {
486: Vec f;
487: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
488: PetscCall(SNESGetFunction(snes, &f, 0, 0));
489: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
490: }
492: PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
494: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
495: Free work space. All PETSc objects should be destroyed when they
496: are no longer needed.
497: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
498: PetscCall(VecDestroy(&x));
499: PetscCall(VecDestroy(&r));
500: PetscCall(MatDestroy(&J));
501: PetscCall(SNESDestroy(&snes));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: static PetscErrorCode FormFunction1_block(SNES snes, Vec x, Vec f, void *dummy)
506: {
507: Vec *xx, *ff, x1, x2, f1, f2;
508: PetscScalar ff_0, ff_1;
509: PetscScalar xx_0, xx_1;
510: PetscInt index, nb;
512: PetscFunctionBeginUser;
513: /* get blocks for function */
514: PetscCall(VecNestGetSubVecs(f, &nb, &ff));
515: f1 = ff[0];
516: f2 = ff[1];
518: /* get blocks for solution */
519: PetscCall(VecNestGetSubVecs(x, &nb, &xx));
520: x1 = xx[0];
521: x2 = xx[1];
523: /* get solution values */
524: index = 0;
525: PetscCall(VecGetValues(x1, 1, &index, &xx_0));
526: PetscCall(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: PetscCall(VecSetValue(f1, index, ff_0, INSERT_VALUES));
535: PetscCall(VecSetValue(f2, index, ff_1, INSERT_VALUES));
537: PetscCall(VecAssemblyBegin(f));
538: PetscCall(VecAssemblyEnd(f));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: static PetscErrorCode FormJacobian1_block(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
543: {
544: Vec *xx, x1, x2;
545: PetscScalar xx_0, xx_1;
546: PetscInt index, nb;
547: PetscScalar A_00, A_01, A_10, A_11;
548: Mat j11, j12, j21, j22;
549: Mat **mats;
551: PetscFunctionBeginUser;
552: /* get blocks for solution */
553: PetscCall(VecNestGetSubVecs(x, &nb, &xx));
554: x1 = xx[0];
555: x2 = xx[1];
557: /* get solution values */
558: index = 0;
559: PetscCall(VecGetValues(x1, 1, &index, &xx_0));
560: PetscCall(VecGetValues(x2, 1, &index, &xx_1));
562: /* get block matrices */
563: PetscCall(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: PetscCall(MatSetValue(j11, 0, 0, A_00, INSERT_VALUES));
577: PetscCall(MatSetValue(j12, 0, 0, A_01, INSERT_VALUES));
578: PetscCall(MatSetValue(j21, 0, 0, A_10, INSERT_VALUES));
579: PetscCall(MatSetValue(j22, 0, 0, A_11, INSERT_VALUES));
581: /* Assemble sub matrix */
582: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
583: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: static PetscErrorCode FormFunction2_block(SNES snes, Vec x, Vec f, void *dummy)
588: {
589: PetscScalar *ff;
590: const PetscScalar *xx;
592: PetscFunctionBeginUser;
593: /*
594: Get pointers to vector data.
595: - For default PETSc vectors, VecGetArray() returns a pointer to
596: the data array. Otherwise, the routine is implementation dependent.
597: - You MUST call VecRestoreArray() when you no longer need access to
598: the array.
599: */
600: PetscCall(VecGetArrayRead(x, &xx));
601: PetscCall(VecGetArray(f, &ff));
603: /*
604: Compute function
605: */
606: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
607: ff[1] = xx[1];
609: /*
610: Restore vectors
611: */
612: PetscCall(VecRestoreArrayRead(x, &xx));
613: PetscCall(VecRestoreArray(f, &ff));
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: static PetscErrorCode FormJacobian2_block(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
618: {
619: const PetscScalar *xx;
620: PetscScalar A[4];
621: PetscInt idx[2] = {0, 1};
623: PetscFunctionBeginUser;
624: /*
625: Get pointer to vector data
626: */
627: PetscCall(VecGetArrayRead(x, &xx));
629: /*
630: Compute Jacobian entries and insert into matrix.
631: - Since this is such a small problem, we set all entries for
632: the matrix at once.
633: */
634: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
635: A[1] = 0.0;
636: A[2] = 0.0;
637: A[3] = 1.0;
638: PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
640: /*
641: Restore vector
642: */
643: PetscCall(VecRestoreArrayRead(x, &xx));
645: /*
646: Assemble matrix
647: */
648: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
649: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: int main(int argc, char **argv)
654: {
655: PetscMPIInt size;
657: PetscFunctionBeginUser;
658: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
659: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
660: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
661: PetscCall(assembled_system());
662: PetscCall(block_system());
663: PetscCall(PetscFinalize());
664: return 0;
665: }
667: /*TEST
669: test:
670: args: -snes_monitor_short
671: requires: !single
673: TEST*/