1 static char help[] = "Basic vector routines.\n\n";
2
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 */
9
10 #include <petscvec.h>
11
main(int argc,char ** argv)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;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
23
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));
33
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));
40
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));
61
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 */
67
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));
71
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"));
75
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));
81
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));
87
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));
93
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));
99
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));
109
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));
115
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));
121
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));
127
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));
136
137 dots[0] = one;
138 dots[1] = three;
139 dots[2] = two;
140
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));
153
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 }
165
166 /*TEST
167
168 testset:
169 output_file: output/ex1_1.out
170 # This is a test where the exact numbers are critical
171 diff_args: -j
172
173 test:
174
175 test:
176 suffix: cuda
177 args: -vec_type cuda
178 requires: cuda
179
180 test:
181 suffix: kokkos
182 args: -vec_type kokkos
183 requires: kokkos_kernels
184
185 test:
186 suffix: hip
187 args: -vec_type hip
188 requires: hip
189
190 test:
191 suffix: 2
192 nsize: 2
193
194 test:
195 suffix: 2_cuda
196 nsize: 2
197 args: -vec_type cuda
198 requires: cuda
199
200 test:
201 suffix: 2_kokkos
202 nsize: 2
203 args: -vec_type kokkos
204 requires: kokkos_kernels
205
206 test:
207 suffix: 2_hip
208 nsize: 2
209 args: -vec_type hip
210 requires: hip
211
212 TEST*/
213