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