Actual source code: ex1.c
1: static char help[] = "Basic vector routines.\n\n";
3: /*
4: Include "petscvec.h" so that we can use vectors. Note that this file
5: automatically includes:
6: petscsys.h - base PETSc routines petscis.h - index sets
7: petscviewer.h - viewers
8: */
10: #include <petscvec.h>
12: int main(int argc, char **argv)
13: {
14: Vec x, y, w; /* vectors */
15: Vec *z; /* array of vectors */
16: PetscReal norm, v, v1, v2, maxval;
17: PetscInt n = 20, maxind;
18: PetscScalar one = 1.0, two = 2.0, three = 3.0, dots[3], dot;
20: PetscFunctionBeginUser;
21: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
22: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
24: /*
25: Create a vector, specifying only its global dimension.
26: When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector format
27: (currently parallel, shared, or sequential) is determined at runtime. Also, the
28: parallel partitioning of the vector is determined by PETSc at runtime.
29: */
30: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
31: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
32: PetscCall(VecSetFromOptions(x));
34: /*
35: Duplicate some work vectors (of the same format and
36: partitioning as the initial vector).
37: */
38: PetscCall(VecDuplicate(x, &y));
39: PetscCall(VecDuplicate(x, &w));
41: /*
42: Duplicate more work vectors (of the same format and
43: partitioning as the initial vector). Here we duplicate
44: an array of vectors, which is often more convenient than
45: duplicating individual ones.
46: */
47: PetscCall(VecDuplicateVecs(x, 3, &z));
48: /*
49: Set the vectors to entries to a constant value.
50: */
51: PetscCall(VecSet(x, one));
52: PetscCall(VecSet(y, two));
53: PetscCall(VecSet(z[0], one));
54: PetscCall(VecSet(z[1], two));
55: PetscCall(VecSet(z[2], three));
56: /*
57: Demonstrate various basic vector routines.
58: */
59: PetscCall(VecDot(x, y, &dot));
60: PetscCall(VecMDot(x, 3, z, dots));
62: /*
63: Note: If using a complex numbers version of PETSc, then
64: PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
65: (when using real numbers) it is undefined.
66: */
68: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Vector length %" PetscInt_FMT "\n", n));
69: PetscCall(VecMax(x, &maxind, &maxval));
70: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecMax %g, VecInd %" PetscInt_FMT "\n", (double)maxval, maxind));
72: PetscCall(VecMin(x, &maxind, &maxval));
73: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecMin %g, VecInd %" PetscInt_FMT "\n", (double)maxval, maxind));
74: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "All other values should be near zero\n"));
76: PetscCall(VecScale(x, two));
77: PetscCall(VecNorm(x, NORM_2, &norm));
78: v = norm - 2.0 * PetscSqrtReal((PetscReal)n);
79: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
80: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecScale %g\n", (double)v));
82: PetscCall(VecCopy(x, w));
83: PetscCall(VecNorm(w, NORM_2, &norm));
84: v = norm - 2.0 * PetscSqrtReal((PetscReal)n);
85: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
86: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecCopy %g\n", (double)v));
88: PetscCall(VecAXPY(y, three, x));
89: PetscCall(VecNorm(y, NORM_2, &norm));
90: v = norm - 8.0 * PetscSqrtReal((PetscReal)n);
91: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
92: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecAXPY %g\n", (double)v));
94: PetscCall(VecAYPX(y, two, x));
95: PetscCall(VecNorm(y, NORM_2, &norm));
96: v = norm - 18.0 * PetscSqrtReal((PetscReal)n);
97: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
98: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecAYPX %g\n", (double)v));
100: PetscCall(VecSwap(x, y));
101: PetscCall(VecNorm(y, NORM_2, &norm));
102: v = norm - 2.0 * PetscSqrtReal((PetscReal)n);
103: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
104: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecSwap %g\n", (double)v));
105: PetscCall(VecNorm(x, NORM_2, &norm));
106: v = norm - 18.0 * PetscSqrtReal((PetscReal)n);
107: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
108: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecSwap %g\n", (double)v));
110: PetscCall(VecWAXPY(w, two, x, y));
111: PetscCall(VecNorm(w, NORM_2, &norm));
112: v = norm - 38.0 * PetscSqrtReal((PetscReal)n);
113: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
114: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecWAXPY %g\n", (double)v));
116: PetscCall(VecPointwiseMult(w, y, x));
117: PetscCall(VecNorm(w, NORM_2, &norm));
118: v = norm - 36.0 * PetscSqrtReal((PetscReal)n);
119: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
120: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecPointwiseMult %g\n", (double)v));
122: PetscCall(VecPointwiseDivide(w, x, y));
123: PetscCall(VecNorm(w, NORM_2, &norm));
124: v = norm - 9.0 * PetscSqrtReal((PetscReal)n);
125: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
126: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecPointwiseDivide %g\n", (double)v));
128: PetscCall(VecSetValue(y, 0, 0.0, INSERT_VALUES));
129: PetscCall(VecAssemblyBegin(y));
130: PetscCall(VecAssemblyEnd(y));
131: PetscCall(VecPointwiseDivide(w, x, y));
132: PetscCall(VecNorm(w, NORM_2, &norm));
133: v = norm - 9.0 * PetscSqrtReal((PetscReal)n);
134: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
135: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecPointwiseDivide %g\n", (double)v));
137: dots[0] = one;
138: dots[1] = three;
139: dots[2] = two;
141: PetscCall(VecSet(x, one));
142: PetscCall(VecMAXPY(x, 3, dots, z));
143: PetscCall(VecNorm(z[0], NORM_2, &norm));
144: v = norm - PetscSqrtReal((PetscReal)n);
145: if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
146: PetscCall(VecNorm(z[1], NORM_2, &norm));
147: v1 = norm - 2.0 * PetscSqrtReal((PetscReal)n);
148: if (v1 > -PETSC_SMALL && v1 < PETSC_SMALL) v1 = 0.0;
149: PetscCall(VecNorm(z[2], NORM_2, &norm));
150: v2 = norm - 3.0 * PetscSqrtReal((PetscReal)n);
151: if (v2 > -PETSC_SMALL && v2 < PETSC_SMALL) v2 = 0.0;
152: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecMAXPY %g %g %g \n", (double)v, (double)v1, (double)v2));
154: /*
155: Free work space. All PETSc objects should be destroyed when they
156: are no longer needed.
157: */
158: PetscCall(VecDestroy(&x));
159: PetscCall(VecDestroy(&y));
160: PetscCall(VecDestroy(&w));
161: PetscCall(VecDestroyVecs(3, &z));
162: PetscCall(PetscFinalize());
163: return 0;
164: }
166: /*TEST
168: testset:
169: output_file: output/ex1_1.out
170: # This is a test where the exact numbers are critical
171: diff_args: -j
173: test:
175: test:
176: suffix: cuda
177: args: -vec_type cuda
178: requires: cuda
180: test:
181: suffix: kokkos
182: args: -vec_type kokkos
183: requires: kokkos_kernels
185: test:
186: suffix: hip
187: args: -vec_type hip
188: requires: hip
190: test:
191: suffix: 2
192: nsize: 2
194: test:
195: suffix: 2_cuda
196: nsize: 2
197: args: -vec_type cuda
198: requires: cuda
200: test:
201: suffix: 2_kokkos
202: nsize: 2
203: args: -vec_type kokkos
204: requires: kokkos_kernels
206: test:
207: suffix: 2_hip
208: nsize: 2
209: args: -vec_type hip
210: requires: hip
212: TEST*/