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;
23: PetscInitialize(&argc, &argv, (char *)0, help);
25: /* create vectors X,Y and F and set values in it*/
26: VecCreate(PETSC_COMM_SELF, &X);
27: VecSetSizes(X, N, N);
28: VecSetFromOptions(X);
29: VecDuplicate(X, &Y);
30: VecDuplicate(X, &F);
31: PetscObjectSetName((PetscObject)F, "F");
32: for (i = 0; i < N; i++) {
33: value = i;
34: VecSetValues(X, 1, &i, &value, INSERT_VALUES);
35: value = 100 + i;
36: VecSetValues(Y, 1, &i, &value, INSERT_VALUES);
37: }
38: VecSet(F, zero);
40: /* Create subvectors X1,X2,Y1,Y2,F1,F2 */
41: VecCreate(PETSC_COMM_SELF, &X1);
42: VecSetSizes(X1, N / 2, N / 2);
43: VecSetFromOptions(X1);
44: VecDuplicate(X1, &X2);
45: VecDuplicate(X1, &Y1);
46: VecDuplicate(X1, &Y2);
47: VecDuplicate(X1, &F1);
48: VecDuplicate(X1, &F2);
50: /* Get array pointers for X,Y,F */
51: VecGetArrayRead(X, &x);
52: VecGetArrayRead(Y, &y);
53: VecGetArray(F, &f);
54: /* Share X,Y,F array pointers with subvectors */
55: VecPlaceArray(X1, x);
56: VecPlaceArray(X2, x + N / 2);
57: VecPlaceArray(Y1, y);
58: VecPlaceArray(Y2, y + N / 2);
59: VecPlaceArray(F1, f);
60: VecPlaceArray(F2, f + N / 2);
62: /* Do subvector addition */
63: VecWAXPY(F1, 1.0, X1, Y1);
64: VecWAXPY(F2, 1.0, X2, Y2);
66: /* Reset subvectors */
67: VecResetArray(X1);
68: VecResetArray(X2);
69: VecResetArray(Y1);
70: VecResetArray(Y2);
71: VecResetArray(F1);
72: VecResetArray(F2);
74: /* Restore X,Y,and F */
75: VecRestoreArrayRead(X, &x);
76: VecRestoreArrayRead(Y, &y);
77: VecRestoreArray(F, &f);
79: PetscPrintf(PETSC_COMM_SELF, "F = X + Y\n");
80: VecView(F, 0);
81: /* Destroy vectors */
82: VecDestroy(&X);
83: VecDestroy(&Y);
84: VecDestroy(&F);
85: VecDestroy(&X1);
86: VecDestroy(&Y1);
87: VecDestroy(&F1);
88: VecDestroy(&X2);
89: VecDestroy(&Y2);
90: VecDestroy(&F2);
92: PetscFinalize();
93: return 0;
94: }