Actual source code: ex39.c
1: static char help[] = "This example is intended for showing how subvectors can\n\
2: share the pointer with the main vector using VecGetArray()\n\
3: and VecPlaceArray() routines so that vector operations done\n\
4: on the subvectors automatically modify the values in the main vector.\n\n";
6: #include <petscvec.h>
8: /* This example shares the array pointers of vectors X,Y,and F with subvectors
9: X1,X2,Y1,Y2,F1,F2 and does vector addition on the subvectors F1 = X1 + Y1, F2 = X2 + Y2 so
10: that F gets updated as a result of sharing the pointers.
11: */
13: int main(int argc, char **argv)
14: {
15: PetscInt N = 10, i;
16: Vec X, Y, F, X1, Y1, X2, Y2, F1, F2;
17: PetscScalar value, zero = 0.0;
18: const PetscScalar *x, *y;
19: PetscScalar *f;
21: PetscFunctionBeginUser;
22: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
24: /* create vectors X,Y and F and set values in it*/
25: PetscCall(VecCreate(PETSC_COMM_SELF, &X));
26: PetscCall(VecSetSizes(X, N, N));
27: PetscCall(VecSetFromOptions(X));
28: PetscCall(VecDuplicate(X, &Y));
29: PetscCall(VecDuplicate(X, &F));
30: PetscCall(PetscObjectSetName((PetscObject)F, "F"));
31: for (i = 0; i < N; i++) {
32: value = i;
33: PetscCall(VecSetValues(X, 1, &i, &value, INSERT_VALUES));
34: value = 100 + i;
35: PetscCall(VecSetValues(Y, 1, &i, &value, INSERT_VALUES));
36: }
37: PetscCall(VecSet(F, zero));
39: /* Create subvectors X1,X2,Y1,Y2,F1,F2 */
40: PetscCall(VecCreate(PETSC_COMM_SELF, &X1));
41: PetscCall(VecSetSizes(X1, N / 2, N / 2));
42: PetscCall(VecSetFromOptions(X1));
43: PetscCall(VecDuplicate(X1, &X2));
44: PetscCall(VecDuplicate(X1, &Y1));
45: PetscCall(VecDuplicate(X1, &Y2));
46: PetscCall(VecDuplicate(X1, &F1));
47: PetscCall(VecDuplicate(X1, &F2));
49: /* Get array pointers for X,Y,F */
50: PetscCall(VecGetArrayRead(X, &x));
51: PetscCall(VecGetArrayRead(Y, &y));
52: PetscCall(VecGetArray(F, &f));
53: /* Share X,Y,F array pointers with subvectors */
54: PetscCall(VecPlaceArray(X1, x));
55: PetscCall(VecPlaceArray(X2, x + N / 2));
56: PetscCall(VecPlaceArray(Y1, y));
57: PetscCall(VecPlaceArray(Y2, y + N / 2));
58: PetscCall(VecPlaceArray(F1, f));
59: PetscCall(VecPlaceArray(F2, f + N / 2));
61: /* Do subvector addition */
62: PetscCall(VecWAXPY(F1, 1.0, X1, Y1));
63: PetscCall(VecWAXPY(F2, 1.0, X2, Y2));
65: /* Reset subvectors */
66: PetscCall(VecResetArray(X1));
67: PetscCall(VecResetArray(X2));
68: PetscCall(VecResetArray(Y1));
69: PetscCall(VecResetArray(Y2));
70: PetscCall(VecResetArray(F1));
71: PetscCall(VecResetArray(F2));
73: /* Restore X,Y,and F */
74: PetscCall(VecRestoreArrayRead(X, &x));
75: PetscCall(VecRestoreArrayRead(Y, &y));
76: PetscCall(VecRestoreArray(F, &f));
78: PetscCall(PetscPrintf(PETSC_COMM_SELF, "F = X + Y\n"));
79: PetscCall(VecView(F, 0));
80: /* Destroy vectors */
81: PetscCall(VecDestroy(&X));
82: PetscCall(VecDestroy(&Y));
83: PetscCall(VecDestroy(&F));
84: PetscCall(VecDestroy(&X1));
85: PetscCall(VecDestroy(&Y1));
86: PetscCall(VecDestroy(&F1));
87: PetscCall(VecDestroy(&X2));
88: PetscCall(VecDestroy(&Y2));
89: PetscCall(VecDestroy(&F2));
91: PetscCall(PetscFinalize());
92: return 0;
93: }