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