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*/