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:   dots[0] = one;
129:   dots[1] = three;
130:   dots[2] = two;

132:   PetscCall(VecSet(x, one));
133:   PetscCall(VecMAXPY(x, 3, dots, z));
134:   PetscCall(VecNorm(z[0], NORM_2, &norm));
135:   v = norm - PetscSqrtReal((PetscReal)n);
136:   if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
137:   PetscCall(VecNorm(z[1], NORM_2, &norm));
138:   v1 = norm - 2.0 * PetscSqrtReal((PetscReal)n);
139:   if (v1 > -PETSC_SMALL && v1 < PETSC_SMALL) v1 = 0.0;
140:   PetscCall(VecNorm(z[2], NORM_2, &norm));
141:   v2 = norm - 3.0 * PetscSqrtReal((PetscReal)n);
142:   if (v2 > -PETSC_SMALL && v2 < PETSC_SMALL) v2 = 0.0;
143:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecMAXPY %g %g %g \n", (double)v, (double)v1, (double)v2));

145:   /*
146:      Free work space.  All PETSc objects should be destroyed when they
147:      are no longer needed.
148:   */
149:   PetscCall(VecDestroy(&x));
150:   PetscCall(VecDestroy(&y));
151:   PetscCall(VecDestroy(&w));
152:   PetscCall(VecDestroyVecs(3, &z));
153:   PetscCall(PetscFinalize());
154:   return 0;
155: }

157: /*TEST

159:   testset:
160:     output_file: output/ex1_1.out
161:     # This is a test where the exact numbers are critical
162:     diff_args: -j

164:     test:

166:     test:
167:         suffix: cuda
168:         args: -vec_type cuda
169:         requires: cuda

171:     test:
172:         suffix: kokkos
173:         args: -vec_type kokkos
174:         requires: kokkos_kernels

176:     test:
177:         suffix: hip
178:         args: -vec_type hip
179:         requires: hip

181:     test:
182:         suffix: 2
183:         nsize: 2

185:     test:
186:         suffix: 2_cuda
187:         nsize: 2
188:         args: -vec_type cuda
189:         requires: cuda

191:     test:
192:         suffix: 2_kokkos
193:         nsize: 2
194:         args: -vec_type kokkos
195:         requires: kokkos_kernels

197:     test:
198:         suffix: 2_hip
199:         nsize: 2
200:         args: -vec_type hip
201:         requires: hip

203: TEST*/