text: ///
  /// Copyright (c) 2019, Intel Corporation
  ///
  /// Redistribution and use in source and binary forms, with or without
  /// modification, are permitted provided that the following conditions
  /// are met:
  ///
  /// * Redistributions of source code must retain the above copyright
  ///       notice, this list of conditions and the following disclaimer.
  /// * Redistributions in binary form must reproduce the above
  ///       copyright notice, this list of conditions and the following
  ///       disclaimer in the documentation and/or other materials provided
  ///       with the distribution.
  /// * Neither the name of Intel Corporation nor the names of its
  ///       contributors may be used to endorse or promote products
  ///       derived from this software without specific prior written
  ///       permission.
  ///
  /// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  /// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  /// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  /// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  /// COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  /// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  /// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  /// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  /// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  /// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  /// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  /// POSSIBILITY OF SUCH DAMAGE.
  
  //////////////////////////////////////////////////////////////////////
  ///
  /// NAME:    nstream
  ///
  /// PURPOSE: To compute memory bandwidth when adding a vector of a given
  ///          number of double precision values to the scalar multiple of
  ///          another vector of the same length, and storing the result in
  ///          a third vector.
  ///
  /// USAGE:   The program takes as input the number
  ///          of iterations to loop over the triad vectors and
  ///          the length of the vectors.
  ///
  ///          <progname> <# iterations> <vector length>
  ///
  ///          The output consists of diagnostics to make sure the
  ///          algorithm worked, and of timing statistics.
  ///
  /// NOTES:   Bandwidth is determined as the number of words read, plus the
  ///          number of words written, times the size of the words, divided
  ///          by the execution time. For a vector length of N, the total
  ///          number of words read and written is 4*N*sizeof(double).
  ///
  ///
  /// HISTORY: This code is loosely based on the Stream benchmark by John
  ///          McCalpin, but does not follow all the Stream rules. Hence,
  ///          reported results should not be associated with Stream in
  ///          external publications
  ///
  ///          Converted to C++11 by Jeff Hammond, November 2017.
  ///          Converted to C11 by Jeff Hammond, February 2019.
  ///
  //////////////////////////////////////////////////////////////////////
  
  #pragma omp requires unified_shared_memory
  
  #include "prk_util.h"
  #include "prk_openmp.h"
  
  int main(int argc, char * argv[])
  {
    printf("Parallel Research Kernels version %d\n", PRKVERSION );
    printf("C11/OpenMP TARGET STREAM triad: A = B + scalar * C\n");
  
    //////////////////////////////////////////////////////////////////////
    /// Read and test input parameters
    //////////////////////////////////////////////////////////////////////
  
    if (argc < 3) {
      printf("Usage: <# iterations> <vector length>\n");
      return 1;
    }
  
    int iterations = atoi(argv[1]);
    if (iterations < 1) {
      printf("ERROR: iterations must be >= 1\n");
      return 1;
    }
  
    // length of a the vector
    size_t length = atol(argv[2]);
    if (length <= 0) {
      printf("ERROR: Vector length must be greater than 0\n");
      return 1;
    }
  
    int device = (argc > 3) ? atol(argv[3]) : omp_get_initial_device();
    if ( (device < 0 || omp_get_num_devices() <= device ) && (device != omp_get_initial_device()) ) {
      printf("ERROR: device number %d is not valid.\n", device);
      return 1;
    }
  
    printf("Number of iterations = %d\n", iterations);
    printf("Vector length        = %zu\n", length);
    printf("OpenMP Device        = %d\n", device);
  
    //////////////////////////////////////////////////////////////////////
    // Allocate space and perform the computation
    //////////////////////////////////////////////////////////////////////
  
    double nstream_time = 0.0;
  
    int host = omp_get_initial_device();
  
    size_t bytes = length*sizeof(double);
    double * restrict A = omp_target_alloc(bytes, host);
    double * restrict B = omp_target_alloc(bytes, host);
    double * restrict C = omp_target_alloc(bytes, host);
  
    double scalar = 3.0;
  
    #pragma omp parallel for simd schedule(static)
    for (size_t i=0; i<length; i++) {
        A[i] = 0.0;
        B[i] = 2.0;
        C[i] = 2.0;
    }
  
    {
      for (int iter = 0; iter<=iterations; iter++) {
  
        if (iter==1) nstream_time = omp_get_wtime();
  
        #pragma omp target teams distribute parallel for simd schedule(static) device(device)
        for (size_t i=0; i<length; i++) {
            A[i] += B[i] + scalar * C[i];
        }
      }
      nstream_time = omp_get_wtime() - nstream_time;
    }
  
    omp_target_free(C, host);
    omp_target_free(B, host);
  
    //////////////////////////////////////////////////////////////////////
    /// Analyze and output results
    //////////////////////////////////////////////////////////////////////
  
    double ar = 0.0;
    double br = 2.0;
    double cr = 2.0;
    for (int i=0; i<=iterations; i++) {
        ar += br + scalar * cr;
    }
  
    ar *= length;
  
    double asum = 0.0;
    #pragma omp parallel for reduction(+:asum)
    for (size_t i=0; i<length; i++) {
        asum += fabs(A[i]);
    }
  
    omp_target_free(A, host);
  
    double epsilon=1.e-8;
    if (fabs(ar-asum)/asum > epsilon) {
        printf("Failed Validation on output array\n"
               "       Expected checksum: %lf\n"
               "       Observed checksum: %lf\n"
               "ERROR: solution did not validate\n", ar, asum);
        return 1;
    } else {
        printf("Solution validates\n");
        double avgtime = nstream_time/iterations;
        double nbytes = 4.0 * length * sizeof(double);
        printf("Rate (MB/s): %lf Avg time (s): %lf\n", 1.e-6*nbytes/avgtime, avgtime);
    }
  
    return 0;
  }
  