#include <cassert>

#include <cmath>

#include <cstdio>

#include <cstdlib>

#include "zfparray2.h"

#include "zfparray3.h"

#define EPSILON 1e-3

// random integer in {begin, ..., end - 1}

static size_t

rand(size_t begin, size_t end)

{

return begin + size_t(rand()) % (end - begin);

}

// ensure f and g are sufficiently close

static void

verify(double f, double g)

{

if (std::fabs(f - g) > EPSILON) {

fprintf(stderr, "error: %g != %g\n", f, g);

exit(EXIT_FAILURE);

}

}

int main(int argc, char* argv[])

{

size_t nx = 16;

size_t ny = 16;

size_t nz = 16;

size_t x0 = rand(0, nx);

size_t y0 = rand(0, ny);

size_t z0 = rand(0, nz);

size_t mx = rand(1, nx - x0);

size_t my = rand(1, ny - y0);

size_t mz = rand(1, nz - z0);

double rate = 16;

// Usage: test [nx ny nz [x0 y0 z0 mx my mz]]

switch (argc) {

case 10:

if (sscanf(argv[4], "%zu", &x0) != 1 ||

sscanf(argv[5], "%zu", &y0) != 1 ||

sscanf(argv[6], "%zu", &z0) != 1 ||

sscanf(argv[7], "%zu", &mx) != 1 ||

sscanf(argv[8], "%zu", &my) != 1 ||

sscanf(argv[9], "%zu", &mz) != 1)

return EXIT_FAILURE;

// FALLTHROUGH

case 4:

if (sscanf(argv[1], "%zu", &nx) != 1 ||

sscanf(argv[2], "%zu", &ny) != 1 ||

sscanf(argv[3], "%zu", &nz) != 1)

return EXIT_FAILURE;

// FALLTHROUGH

case 1:

break;

}

printf("a(%zu, %zu, %zu)\n", nx, ny, nz);

printf("v(%zu, %zu, %zu) + (%zu, %zu, %zu)\n", mx, my, mz, x0, y0, z0);

// initialize 3D array to linear function

zfp::array3<double> a(nx, ny, nz, rate);

for (size_t z = 0; z < nz; z++)

for (size_t y = 0; y < ny; y++)

for (size_t x = 0; x < nx; x++)

a(x, y, z) = static_cast<double>(x + nx * (y + ny * z));

// rectangular view into a

printf("\n3D view\n");

zfp::array3<double>::view v(&a, x0, y0, z0, mx, my, mz);

for (size_t z = 0; z < v.size_z(); z++)

for (size_t y = 0; y < v.size_y(); y++)

for (size_t x = 0; x < v.size_x(); x++) {

printf("%zu %zu %zu: %g %g\n", x, y, z, (double)a(x0 + x, y0 + y, z0 + z), (double)v(x, y, z));

verify(a(x0 + x, y0 + y, z0 + z), v(x, y, z));

}

// flat view of all of a

printf("\n3D flat view\n");

zfp::array3<double>::flat_view fv(&a);

for (size_t z = 0; z < fv.size_z(); z++)

for (size_t y = 0; y < fv.size_y(); y++)

for (size_t x = 0; x < fv.size_x(); x++) {

printf("%zu %zu %zu: %g %g\n", x, y, z, (double)a(x, y, z), (double)fv[fv.index(x, y, z)]);

verify(a(x, y, z), fv[fv.index(x, y, z)]);

}

// nested view of all of a

printf("\n3D nested view\n");

zfp::array3<double>::nested_view nv(&a);

for (size_t z = 0; z < v.size_z(); z++)

for (size_t y = 0; y < v.size_y(); y++)

for (size_t x = 0; x < v.size_x(); x++) {

printf("%zu %zu %zu: %g %g\n", x, y, z, (double)a(x, y, z), (double)nv[z][y][x]);

verify(a(x, y, z), nv[z][y][x]);

}

// pointers and iterators into a via view v

printf("\n3D view pointers and iterators\n");

zfp::array3<double>::view::const_reference vr = v(0, 0, 0);

zfp::array3<double>::view::const_pointer p = &vr;

p = &v(0, 0, 0);

for (zfp::array3<double>::view::const_iterator it = v.begin(); it != v.end(); it++) {

size_t x = it.i();

size_t y = it.j();

size_t z = it.k();

verify(*it, p[x + mx * (y + my * z)]);

}

// pointers and iterators into a via flat view fv

printf("\n3D flat view pointers and iterators\n");

zfp::array3<double>::flat_view::const_reference fvr = fv[0];

zfp::array3<double>::flat_view::const_pointer fp = &fvr;

fp = &fv(0, 0, 0);

for (zfp::array3<double>::flat_view::const_iterator it = fv.begin(); it != fv.end(); it++) {

size_t x = it.i();

size_t y = it.j();

size_t z = it.k();

verify(*it, fp[x + nx * (y + ny * z)]);

}

// 2D slice of a

printf("\n2D slice\n");

size_t z = rand(0, nv.size_z());

zfp::array3<double>::nested_view2 slice2(nv[z]);

for (size_t y = 0; y < slice2.size_y(); y++)

for (size_t x = 0; x < slice2.size_x(); x++) {

printf("%zu %zu %zu: %g %g\n", x, y, z, (double)a(x, y, z), (double)slice2[y][x]);

verify(a(x, y, z), slice2[y][x]);

}

// 2D array constructed from 2D slice (exercises deep copy via iterator)

printf("\n2D array from 2D slice\n");

zfp::array2<double> b(slice2);

for (size_t y = 0; y < b.size_y(); y++)

for (size_t x = 0; x < b.size_x(); x++) {

printf("%zu %zu: %g %g\n", x, y, (double)b(x, y), (double)slice2[y][x]);

verify(b(x, y), slice2[y][x]);

}

// 1D slice of a

printf("\n1D slice\n");

size_t y = rand(0, slice2.size_y());

zfp::array3<double>::nested_view1 slice1 = slice2[y];

for (size_t x = 0; x < slice1.size_x(); x++) {

printf("%zu %zu %zu: %g %g\n", x, y, z, (double)a(x, y, z), (double)slice1[x]);

verify(a(x, y, z), slice1[x]);

}

// 2D array constructed from 2D slice of 3D array (exercises deep copy via iterator)

printf("\n2D array from 2D slice of 3D array\n");

zfp::array2<double> c(slice2);

for (size_t y = 0; y < c.size_y(); y++)

for (size_t x = 0; x < c.size_x(); x++) {

printf("%zu %zu: %g %g\n", x, y, (double)c(x, y), (double)slice2[y][x]);

verify(c(x, y), slice2[y][x]);

}

// 2D thread-safe view of c

printf("\n2D private view\n");

zfp::array2<double>::private_const_view d(&c);

for (size_t y = 0; y < c.size_y(); y++)

for (size_t x = 0; x < c.size_x(); x++) {

printf("%zu %zu: %g %g\n", x, y, (double)c(x, y), (double)d(x, y));

verify(c(x, y), d(x, y));

}

printf("\nall tests passed\n");

return 0;

}