Actual source code: ex53.c

  1: static const char help[] = "Tests VecShift()\n\n";

  3: #include <petscvec.h>

  5: static PetscBool PetscIsCloseAtTolScalar(PetscScalar l, PetscScalar r, PetscReal atol, PetscReal rtol)
  6: {
  7:   return (PetscBool)(PetscIsCloseAtTol(PetscRealPart(l), PetscRealPart(r), atol, rtol) && PetscIsCloseAtTol(PetscImaginaryPart(l), PetscImaginaryPart(r), atol, rtol));
  8: }

 10: static PetscErrorCode CheckVecShift(Vec v, PetscInt n, PetscScalar *array_copy, PetscScalar shift)
 11: {
 12:   const PetscScalar *array;

 14:   PetscFunctionBegin;
 15:   for (PetscInt i = 0; i < n; ++i) array_copy[i] += shift;
 16:   PetscCall(VecShift(v, shift));
 17:   PetscCall(VecGetArrayRead(v, &array));
 18:   for (PetscInt i = 0; i < n; ++i) {
 19:     const PetscScalar actual = array[i], expected = array_copy[i];

 21:     PetscCheck(PetscIsCloseAtTolScalar(actual, expected, 1e-12, 0.0), PETSC_COMM_SELF, PETSC_ERR_PLIB, "VecShift() returned array[%" PetscInt_FMT "] %g + %gi != expected_array[%" PetscInt_FMT "] %g + %gi", i, (double)PetscRealPart(actual), (double)PetscImaginaryPart(actual), i, (double)PetscRealPart(expected), (double)PetscImaginaryPart(expected));
 22:   }
 23:   PetscCall(VecRestoreArrayRead(v, &array));
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

 27: int main(int argc, char **argv)
 28: {
 29:   Vec                x;
 30:   PetscInt           n;
 31:   const PetscScalar *array;
 32:   PetscScalar       *array_copy;
 33:   PetscReal          norm_before, norm_after;
 34:   PetscBool          available;

 36:   PetscFunctionBeginUser;
 37:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 39:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 40:   PetscCall(VecSetSizes(x, PETSC_DECIDE, 10));
 41:   PetscCall(VecSetFromOptions(x));

 43:   PetscCall(VecZeroEntries(x));

 45:   // get a copy of the vectors array, anything we do to the vector via VecShift we will also do
 46:   // to the copy, and hence they should always match
 47:   PetscCall(VecGetLocalSize(x, &n));
 48:   PetscCall(PetscMalloc1(n, &array_copy));
 49:   PetscCall(VecGetArrayRead(x, &array));
 50:   PetscCall(PetscArraycpy(array_copy, array, n));
 51:   PetscCall(VecRestoreArrayRead(x, &array));

 53:   PetscCall(CheckVecShift(x, n, array_copy, 0.0));
 54:   PetscCall(CheckVecShift(x, n, array_copy, 1.0));
 55:   PetscCall(CheckVecShift(x, n, array_copy, -1.0));
 56:   PetscCall(CheckVecShift(x, n, array_copy, 15.0));

 58:   PetscCall(VecNorm(x, NORM_2, &norm_before));
 59:   PetscCall(VecNormAvailable(x, NORM_2, &available, &norm_after));
 60:   PetscCheck(available, PETSC_COMM_SELF, PETSC_ERR_PLIB, "VecNormAvailable() returned FALSE right after calling VecNorm()");
 61:   // a shift of zero should not invalidate norms
 62:   PetscCall(CheckVecShift(x, n, array_copy, 0.0));
 63:   PetscCall(VecNormAvailable(x, NORM_2, &available, &norm_after));
 64:   PetscCheck(available, PETSC_COMM_SELF, PETSC_ERR_PLIB, "VecNormAvailable() returned FALSE after calling VecShift() with a shift of 0.0!");
 65:   // these can be compared with equality as the number should not change *at all*
 66:   PetscCheck(norm_before == norm_after, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Norms differ before and after calling VecShift() with shift of 0.0: before %g after %g", (double)norm_before, (double)norm_after);

 68:   PetscCall(PetscFree(array_copy));
 69:   PetscCall(VecDestroy(&x));
 70:   PetscCall(PetscFinalize());
 71:   return 0;
 72: }

 74: /*TEST

 76:   testset:
 77:     output_file: ./output/empty.out
 78:     nsize: {{1 2}}
 79:     test:
 80:       suffix: standard
 81:     test:
 82:       requires: defined(PETSC_USE_SHARED_MEMORY)
 83:       args: -vec_type shared
 84:       suffix: shared
 85:     test:
 86:       requires: viennacl
 87:       args: -vec_type viennacl
 88:       suffix: viennacl
 89:     test:
 90:       requires: kokkos_kernels
 91:       args: -vec_type kokkos
 92:       suffix: kokkos
 93:     test:
 94:       requires: cuda
 95:       args: -vec_type cuda
 96:       suffix: cuda
 97:     test:
 98:       requires: hip
 99:       args: -vec_type hip
100:       suffix: hip

102: TEST*/