diff options
Diffstat (limited to 'Carpet/CarpetLib/src/restrict_4d_rf2.cc')
-rw-r--r-- | Carpet/CarpetLib/src/restrict_4d_rf2.cc | 141 |
1 files changed, 141 insertions, 0 deletions
diff --git a/Carpet/CarpetLib/src/restrict_4d_rf2.cc b/Carpet/CarpetLib/src/restrict_4d_rf2.cc new file mode 100644 index 000000000..77bf2d28b --- /dev/null +++ b/Carpet/CarpetLib/src/restrict_4d_rf2.cc @@ -0,0 +1,141 @@ +#include <algorithm> +#include <cassert> +#include <cmath> +#include <cstdlib> +#include <iostream> + +#include <cctk.h> +#include <cctk_Parameters.h> + +#include "operator_prototypes_4d.hh" +#include "typeprops.hh" + +using namespace std; + + + +namespace CarpetLib { + + + +#define SRCIND4(i,j,k,l) \ + index4 (srcioff + (i), srcjoff + (j), srckoff + (k), srcloff + (l), \ + srciext, srcjext, srckext, srclext) +#define DSTIND4(i,j,k,l) \ + index4 (dstioff + (i), dstjoff + (j), dstkoff + (k), dstloff + (l), \ + dstiext, dstjext, dstkext, dstlext) + + + + template <typename T> + void + restrict_4d_rf2 (T const * restrict const src, + ivect4 const & restrict srcext, + T * restrict const dst, + ivect4 const & restrict dstext, + ibbox4 const & restrict srcbbox, + ibbox4 const & restrict dstbbox, + ibbox4 const & restrict regbbox) + { + if (any (srcbbox.stride() >= regbbox.stride() or + dstbbox.stride() != regbbox.stride())) + { + CCTK_WARN (0, "Internal error: strides disagree"); + } + + if (any (reffact2 * srcbbox.stride() != dstbbox.stride())) { + CCTK_WARN (0, "Internal error: destination strides are not twice the source strides"); + } + + // This could be handled, but is likely to point to an error + // elsewhere + if (regbbox.empty()) { + CCTK_WARN (0, "Internal error: region extent is empty"); + } + + if (not regbbox.is_contained_in(srcbbox) or + not regbbox.is_contained_in(dstbbox)) + { + cerr << "srcbbox: " << srcbbox << endl + << "dstbbox: " << dstbbox << endl + << "regbbox: " << regbbox << endl; + CCTK_WARN (0, "Internal error: region extent is not contained in array extent"); + } + + if (any (srcext != srcbbox.shape() / srcbbox.stride() or + dstext != dstbbox.shape() / dstbbox.stride())) + { + CCTK_WARN (0, "Internal error: array sizes don't agree with bounding boxes"); + } + + + + ivect4 const regext = regbbox.shape() / regbbox.stride(); + assert (all ((regbbox.lower() - srcbbox.lower()) % srcbbox.stride() == 0)); + ivect4 const srcoff = (regbbox.lower() - srcbbox.lower()) / srcbbox.stride(); + assert (all ((regbbox.lower() - dstbbox.lower()) % dstbbox.stride() == 0)); + ivect4 const dstoff = (regbbox.lower() - dstbbox.lower()) / dstbbox.stride(); + + + + ptrdiff_t const srciext = srcext[0]; + ptrdiff_t const srcjext = srcext[1]; + ptrdiff_t const srckext = srcext[2]; + ptrdiff_t const srclext = srcext[3]; + + ptrdiff_t const dstiext = dstext[0]; + ptrdiff_t const dstjext = dstext[1]; + ptrdiff_t const dstkext = dstext[2]; + ptrdiff_t const dstlext = dstext[3]; + + ptrdiff_t const regiext = regext[0]; + ptrdiff_t const regjext = regext[1]; + ptrdiff_t const regkext = regext[2]; + ptrdiff_t const reglext = regext[3]; + + ptrdiff_t const srcioff = srcoff[0]; + ptrdiff_t const srcjoff = srcoff[1]; + ptrdiff_t const srckoff = srcoff[2]; + ptrdiff_t const srcloff = srcoff[3]; + + ptrdiff_t const dstioff = dstoff[0]; + ptrdiff_t const dstjoff = dstoff[1]; + ptrdiff_t const dstkoff = dstoff[2]; + ptrdiff_t const dstloff = dstoff[3]; + + + + // Loop over coarse region +#pragma omp parallel for + for (int l=0; l<reglext; ++l) { + for (int k=0; k<regkext; ++k) { + for (int j=0; j<regjext; ++j) { + for (int i=0; i<regiext; ++i) { + + dst [DSTIND4(i, j, k, l)] = src [SRCIND4(2*i, 2*j, 2*k, 2*l)]; + + } + } + } + } + + } + + + +#define INSTANTIATE(T) \ + template \ + void \ + restrict_4d_rf2 (T const * restrict const src, \ + ivect4 const & restrict srcext, \ + T * restrict const dst, \ + ivect4 const & restrict dstext, \ + ibbox4 const & restrict srcbbox, \ + ibbox4 const & restrict dstbbox, \ + ibbox4 const & restrict regbbox); +#include "instantiate" +#undef INSTANTIATE + + + +} // namespace CarpetLib |