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