Commit d73de808 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

stage 1

parent 0d8a2890
......@@ -621,7 +621,7 @@ template<typename T, typename T0> aligned_array<T> alloc_tmp
}
template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
const cfmav<Cmplx<T>> &src, Cmplx<native_simd<T>> *MRUTIL_RESTRICT dst)
const fmav<Cmplx<T>> &src, Cmplx<native_simd<T>> *MRUTIL_RESTRICT dst)
{
for (size_t i=0; i<it.length_in(); ++i)
for (size_t j=0; j<vlen; ++j)
......@@ -632,7 +632,7 @@ template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
}
template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
const cfmav<T> &src, native_simd<T> *MRUTIL_RESTRICT dst)
const fmav<T> &src, native_simd<T> *MRUTIL_RESTRICT dst)
{
for (size_t i=0; i<it.length_in(); ++i)
for (size_t j=0; j<vlen; ++j)
......@@ -640,7 +640,7 @@ template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
}
template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
const cfmav<T> &src, T *MRUTIL_RESTRICT dst)
const fmav<T> &src, T *MRUTIL_RESTRICT dst)
{
if (dst == &src[it.iofs(0)]) return; // in-place
for (size_t i=0; i<it.length_in(); ++i)
......@@ -648,7 +648,7 @@ template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
}
template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
const Cmplx<native_simd<T>> *MRUTIL_RESTRICT src, const fmav<Cmplx<T>> &dst)
const Cmplx<native_simd<T>> *MRUTIL_RESTRICT src, fmav<Cmplx<T>> &dst)
{
for (size_t i=0; i<it.length_out(); ++i)
for (size_t j=0; j<vlen; ++j)
......@@ -656,7 +656,7 @@ template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
}
template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
const native_simd<T> *MRUTIL_RESTRICT src, const fmav<T> &dst)
const native_simd<T> *MRUTIL_RESTRICT src, fmav<T> &dst)
{
for (size_t i=0; i<it.length_out(); ++i)
for (size_t j=0; j<vlen; ++j)
......@@ -664,7 +664,7 @@ template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
}
template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
const T *MRUTIL_RESTRICT src, const fmav<T> &dst)
const T *MRUTIL_RESTRICT src, fmav<T> &dst)
{
if (src == &dst[it.oofs(0)]) return; // in-place
for (size_t i=0; i<it.length_out(); ++i)
......@@ -677,7 +677,7 @@ template <typename T> struct add_vec<Cmplx<T>>
template <typename T> using add_vec_t = typename add_vec<T>::type;
template<typename Tplan, typename T, typename T0, typename Exec>
MRUTIL_NOINLINE void general_nd(const cfmav<T> &in, const fmav<T> &out,
MRUTIL_NOINLINE void general_nd(const fmav<T> &in, fmav<T> &out,
const shape_t &axes, T0 fct, size_t nthreads, const Exec & exec,
const bool allow_inplace=true)
{
......@@ -721,9 +721,9 @@ struct ExecC2C
{
bool forward;
template <typename T0, typename T, size_t vlen> void operator () (
const multi_iter<vlen> &it, const cfmav<Cmplx<T0>> &in,
const fmav<Cmplx<T0>> &out, T * buf, const pocketfft_c<T0> &plan, T0 fct) const
template <typename T0, typename T, size_t vlen> void operator() (
const multi_iter<vlen> &it, const fmav<Cmplx<T0>> &in,
fmav<Cmplx<T0>> &out, T *buf, const pocketfft_c<T0> &plan, T0 fct) const
{
copy_input(it, in, buf);
plan.exec(buf, fct, forward);
......@@ -732,7 +732,7 @@ struct ExecC2C
};
template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
const native_simd<T> *MRUTIL_RESTRICT src, const fmav<T> &dst)
const native_simd<T> *MRUTIL_RESTRICT src, fmav<T> &dst)
{
for (size_t j=0; j<vlen; ++j)
dst[it.oofs(j,0)] = src[0][j];
......@@ -740,8 +740,8 @@ template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
for (i=1; i<it.length_out()-1; i+=2, ++i1, --i2)
for (size_t j=0; j<vlen; ++j)
{
dst[it.oofs(j,i1)] = src[i][j]+src[i+1][j];
dst[it.oofs(j,i2)] = src[i][j]-src[i+1][j];
dst[it.oofs(j,i1)] = src[i][j]+src[i+1][j];
dst[it.oofs(j,i2)] = src[i][j]-src[i+1][j];
}
if (i<it.length_out())
for (size_t j=0; j<vlen; ++j)
......@@ -749,7 +749,7 @@ template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
}
template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
const T *MRUTIL_RESTRICT src, const fmav<T> &dst)
const T *MRUTIL_RESTRICT src, fmav<T> &dst)
{
dst[it.oofs(0)] = src[0];
size_t i=1, i1=1, i2=it.length_out()-1;
......@@ -765,7 +765,7 @@ template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
struct ExecHartley
{
template <typename T0, typename T, size_t vlen> void operator () (
const multi_iter<vlen> &it, const cfmav<T0> &in, const fmav<T0> &out,
const multi_iter<vlen> &it, const fmav<T0> &in, fmav<T0> &out,
T * buf, const pocketfft_r<T0> &plan, T0 fct) const
{
copy_input(it, in, buf);
......@@ -781,8 +781,8 @@ struct ExecDcst
bool cosine;
template <typename T0, typename T, typename Tplan, size_t vlen>
void operator () (const multi_iter<vlen> &it, const cfmav<T0> &in,
const fmav <T0> &out, T * buf, const Tplan &plan, T0 fct) const
void operator () (const multi_iter<vlen> &it, const fmav<T0> &in,
fmav <T0> &out, T * buf, const Tplan &plan, T0 fct) const
{
copy_input(it, in, buf);
plan.exec(buf, fct, ortho, type, cosine);
......@@ -791,7 +791,7 @@ struct ExecDcst
};
template<typename T> MRUTIL_NOINLINE void general_r2c(
const cfmav<T> &in, const fmav<Cmplx<T>> &out, size_t axis, bool forward, T fct,
const fmav<T> &in, fmav<Cmplx<T>> &out, size_t axis, bool forward, T fct,
size_t nthreads)
{
auto plan = get_plan<pocketfft_r<T>>(in.shape(axis));
......@@ -846,7 +846,7 @@ template<typename T> MRUTIL_NOINLINE void general_r2c(
}); // end of parallel region
}
template<typename T> MRUTIL_NOINLINE void general_c2r(
const cfmav<Cmplx<T>> &in, const fmav<T> &out, size_t axis, bool forward, T fct,
const fmav<Cmplx<T>> &in, fmav<T> &out, size_t axis, bool forward, T fct,
size_t nthreads)
{
auto plan = get_plan<pocketfft_r<T>>(out.shape(axis));
......@@ -922,7 +922,7 @@ struct ExecR2R
bool r2c, forward;
template <typename T0, typename T, size_t vlen> void operator () (
const multi_iter<vlen> &it, const cfmav<T0> &in, const fmav<T0> &out, T * buf,
const multi_iter<vlen> &it, const fmav<T0> &in, fmav<T0> &out, T * buf,
const pocketfft_r<T0> &plan, T0 fct) const
{
copy_input(it, in, buf);
......@@ -937,18 +937,18 @@ struct ExecR2R
}
};
template<typename T> void c2c(const cfmav<std::complex<T>> &in,
const fmav<std::complex<T>> &out, const shape_t &axes, bool forward,
template<typename T> void c2c(const fmav<std::complex<T>> &in,
fmav<std::complex<T>> &out, const shape_t &axes, bool forward,
T fct, size_t nthreads=1)
{
util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
if (in.size()==0) return;
cfmav<Cmplx<T>> in2(reinterpret_cast<const Cmplx<T> *>(in.data()), in);
fmav<Cmplx<T>> out2(reinterpret_cast<Cmplx<T> *>(out.data()), out);
fmav<Cmplx<T>> in2(reinterpret_cast<const Cmplx<T> *>(in.data()), in);
fmav<Cmplx<T>> out2(reinterpret_cast<Cmplx<T> *>(out.data()), out, out.writable());
general_nd<pocketfft_c<T>>(in2, out2, axes, fct, nthreads, ExecC2C{forward});
}
template<typename T> void dct(const cfmav<T> &in, const fmav<T> &out,
template<typename T> void dct(const fmav<T> &in, fmav<T> &out,
const shape_t &axes, int type, T fct, bool ortho, size_t nthreads=1)
{
if ((type<1) || (type>4)) throw std::invalid_argument("invalid DCT type");
......@@ -963,7 +963,7 @@ template<typename T> void dct(const cfmav<T> &in, const fmav<T> &out,
general_nd<T_dcst23<T>>(in, out, axes, fct, nthreads, exec);
}
template<typename T> void dst(const cfmav<T> &in, const fmav<T> &out,
template<typename T> void dst(const fmav<T> &in, fmav<T> &out,
const shape_t &axes, int type, T fct, bool ortho, size_t nthreads=1)
{
if ((type<1) || (type>4)) throw std::invalid_argument("invalid DST type");
......@@ -977,18 +977,18 @@ template<typename T> void dst(const cfmav<T> &in, const fmav<T> &out,
general_nd<T_dcst23<T>>(in, out, axes, fct, nthreads, exec);
}
template<typename T> void r2c(const cfmav<T> &in,
const fmav<std::complex<T>> &out, size_t axis, bool forward, T fct,
template<typename T> void r2c(const fmav<T> &in,
fmav<std::complex<T>> &out, size_t axis, bool forward, T fct,
size_t nthreads=1)
{
util::sanity_check_cr(out, in, axis);
if (in.size()==0) return;
fmav<Cmplx<T>> out2(reinterpret_cast<Cmplx<T> *>(out.data()), out);
fmav<Cmplx<T>> out2(reinterpret_cast<Cmplx<T> *>(out.data()), out, out.writable());
general_r2c(in, out2, axis, forward, fct, nthreads);
}
template<typename T> void r2c(const cfmav<T> &in,
const fmav<std::complex<T>> &out, const shape_t &axes,
template<typename T> void r2c(const fmav<T> &in,
fmav<std::complex<T>> &out, const shape_t &axes,
bool forward, T fct, size_t nthreads=1)
{
util::sanity_check_cr(out, in, axes);
......@@ -1000,17 +1000,17 @@ template<typename T> void r2c(const cfmav<T> &in,
c2c(out, out, newaxes, forward, T(1), nthreads);
}
template<typename T> void c2r(const cfmav<std::complex<T>> &in,
const fmav<T> &out, size_t axis, bool forward, T fct, size_t nthreads=1)
template<typename T> void c2r(const fmav<std::complex<T>> &in,
fmav<T> &out, size_t axis, bool forward, T fct, size_t nthreads=1)
{
util::sanity_check_cr(in, out, axis);
if (in.size()==0) return;
cfmav<Cmplx<T>> in2(reinterpret_cast<const Cmplx<T> *>(in.data()), in);
fmav<Cmplx<T>> in2(reinterpret_cast<const Cmplx<T> *>(in.data()), in);
general_c2r(in2, out, axis, forward, fct, nthreads);
}
template<typename T> void c2r(const cfmav<std::complex<T>> &in,
const fmav<T> &out, const shape_t &axes, bool forward, T fct,
template<typename T> void c2r(const fmav<std::complex<T>> &in,
fmav<T> &out, const shape_t &axes, bool forward, T fct,
size_t nthreads=1)
{
if (axes.size()==1)
......@@ -1023,8 +1023,8 @@ template<typename T> void c2r(const cfmav<std::complex<T>> &in,
c2r(atmp, out, axes.back(), forward, fct, nthreads);
}
template<typename T> void r2r_fftpack(const cfmav<T> &in,
const fmav<T> &out, const shape_t &axes, bool real2hermitian, bool forward,
template<typename T> void r2r_fftpack(const fmav<T> &in,
fmav<T> &out, const shape_t &axes, bool real2hermitian, bool forward,
T fct, size_t nthreads=1)
{
util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
......@@ -1033,8 +1033,8 @@ template<typename T> void r2r_fftpack(const cfmav<T> &in,
ExecR2R{real2hermitian, forward});
}
template<typename T> void r2r_separable_hartley(const cfmav<T> &in,
const fmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
template<typename T> void r2r_separable_hartley(const fmav<T> &in,
fmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
{
util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
if (in.size()==0) return;
......@@ -1042,8 +1042,8 @@ template<typename T> void r2r_separable_hartley(const cfmav<T> &in,
false);
}
template<typename T> void r2r_genuine_hartley(const cfmav<T> &in,
const fmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
template<typename T> void r2r_genuine_hartley(const fmav<T> &in,
fmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
{
if (axes.size()==1)
return r2r_separable_hartley(in, out, axes, fct, nthreads);
......
......@@ -16,7 +16,7 @@
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* Copyright (C) 2019 Max-Planck-Society
/* Copyright (C) 2019-2020 Max-Planck-Society
Author: Martin Reinecke */
#ifndef MRUTIL_MAV_H
......@@ -42,6 +42,7 @@ class fmav_info
protected:
shape_t shp;
stride_t str;
size_t sz;
static size_t prod(const shape_t &shape)
{
......@@ -53,18 +54,22 @@ class fmav_info
public:
fmav_info(const shape_t &shape_, const stride_t &stride_)
: shp(shape_), str(stride_)
{ MR_assert(shp.size()==str.size(), "dimensions mismatch"); }
: shp(shape_), str(stride_), sz(prod(shp))
{
MR_assert(shp.size()>0, "at least 1D required");
MR_assert(shp.size()==str.size(), "dimensions mismatch");
}
fmav_info(const shape_t &shape_)
: shp(shape_), str(shape_.size())
: shp(shape_), str(shape_.size()), sz(prod(shp))
{
auto ndim = shp.size();
MR_assert(ndim>0, "at least 1D required");
str[ndim-1]=1;
for (size_t i=2; i<=ndim; ++i)
str[ndim-i] = str[ndim-i+1]*ptrdiff_t(shp[ndim-i+1]);
}
size_t ndim() const { return shp.size(); }
size_t size() const { return prod(shp); }
size_t size() const { return sz; }
const shape_t &shape() const { return shp; }
size_t shape(size_t i) const { return shp[i]; }
const stride_t &stride() const { return str; }
......@@ -86,200 +91,197 @@ class fmav_info
{ return shp==other.shp; }
};
template<typename T, size_t ndim> class cmav;
// "mav" stands for "multidimensional array view"
template<typename T> class cfmav: public fmav_info
template<typename T> class membuf
{
protected:
using Tsp = shared_ptr<vector<T>>;
Tsp ptr;
T *d;
const T *d;
bool rw;
public:
cfmav(const T *d_, const shape_t &shp_, const stride_t &str_)
: fmav_info(shp_, str_), d(const_cast<T *>(d_)) {}
cfmav(const T *d_, const shape_t &shp_)
: fmav_info(shp_), d(const_cast<T *>(d_)) {}
cfmav(const shape_t &shp_)
: fmav_info(shp_), ptr(make_unique<vector<T>>(size())),
d(const_cast<T *>(ptr->data())) {}
cfmav(const T* d_, const fmav_info &info)
: fmav_info(info), d(const_cast<T *>(d_)) {}
cfmav(const cfmav &other) = default;
cfmav(cfmav &&other) = default;
membuf(T *d_, bool rw_=false)
: d(d_), rw(rw_) {}
membuf(const T *d_)
: d(d_), rw(false) {}
membuf(size_t sz)
: ptr(make_unique<vector<T>>(sz)), d(ptr->data()), rw(true) {}
membuf(const membuf &other) = default;
membuf(membuf &&other) = default;
// Not for public use!
cfmav(const T *d_, const Tsp &p, const shape_t &shp_, const stride_t &str_)
: fmav_info(shp_, str_), ptr(p), d(const_cast<T *>(d_)) {}
membuf(T *d_, const Tsp &p, bool rw_)
: ptr(p), d(d_), rw(rw_) {}
template<typename I> const T &val(I i) const
{ return d[i]; }
template<typename I> T &val(I i)
{
MR_assert(rw, "array is nor writable");
return const_cast<T *>(d)[i];
}
template<typename I> const T &operator[](I i) const
{ return d[i]; }
template<typename I> T &operator[](I i)
{
MR_assert(rw, "array is nor writable");
return const_cast<T *>(d)[i];
}
const T *data() const
{ return d; }
T *data()
{
MR_assert(rw, "array is nor writable");
return const_cast<T *>(d);
}
bool writable() const { return rw; }
};
template<typename T> class fmav: public cfmav<T>
// "mav" stands for "multidimensional array view"
template<typename T> class fmav: public fmav_info, public membuf<T>
{
protected:
using Tsp = shared_ptr<vector<T>>;
using parent = cfmav<T>;
using parent::d;
using parent::shp;
using parent::str;
using typename membuf<T>::Tsp;
public:
fmav(T *d_, const shape_t &shp_, const stride_t &str_)
: parent(d_, shp_, str_) {}
fmav(T *d_, const shape_t &shp_)
: parent(d_, shp_) {}
fmav(const T *d_, const shape_t &shp_, const stride_t &str_)
: fmav_info(shp_, str_), membuf<T>(d_) {}
fmav(const T *d_, const shape_t &shp_)
: fmav_info(shp_), membuf<T>(d_) {}
fmav(T *d_, const shape_t &shp_, const stride_t &str_, bool rw_)
: fmav_info(shp_, str_), membuf<T>(d_,rw_) {}
fmav(T *d_, const shape_t &shp_, bool rw_)
: fmav_info(shp_), membuf<T>(d_,rw_) {}
fmav(const shape_t &shp_)
: parent(shp_) {}
fmav(T* d_, const fmav_info &info)
: parent(d_, info) {}
fmav(const fmav &other) = default;
: fmav_info(shp_), membuf<T>(size()) {}
fmav(T* d_, const fmav_info &info, bool rw_=false)
: fmav_info(info), membuf<T>(d_, rw_) {}
fmav(const T* d_, const fmav_info &info)
: fmav_info(info), membuf<T>(d_) {}
fmav(const fmav &other) = delete;
fmav(fmav &&other) = default;
// Not for public use!
fmav(T *d_, const Tsp &p, const shape_t &shp_, const stride_t &str_)
: parent(d_, p, shp_, str_) {}
template<typename I> T &operator[](I i) const
{ return d[i]; }
using parent::shape;
using parent::stride;
T *data() const
{ return d; }
using parent::last_contiguous;
using parent::contiguous;
using parent::conformable;
fmav(T *d_, const Tsp &p, const shape_t &shp_, const stride_t &str_, bool rw_)
: fmav_info(shp_, str_), membuf<T>(d_, p, rw_) {}
};
template<typename T, size_t ndim> class cmav
template<size_t ndim> class mav_info
{
static_assert((ndim>0) && (ndim<4), "only supports 1D, 2D, and 3D arrays");
protected:
using Tsp = shared_ptr<vector<T>>;
array<size_t, ndim> shp;
array<ptrdiff_t, ndim> str;
Tsp ptr;
T *d;
using shape_t = array<size_t, ndim>;
using stride_t = array<ptrdiff_t, ndim>;
public:
cmav(const T *d_, const array<size_t,ndim> &shp_,
const array<ptrdiff_t,ndim> &str_)
: shp(shp_), str(str_), d(const_cast<T *>(d_)) {}
cmav(const T *d_, const array<size_t,ndim> &shp_)
: shp(shp_), d(const_cast<T *>(d_))
shape_t shp;
stride_t str;
size_t sz;
static size_t prod(const shape_t &shape)
{
str[ndim-1]=1;
for (size_t i=2; i<=ndim; ++i)
str[ndim-i] = str[ndim-i+1]*shp[ndim-i+1];
size_t res=1;
for (auto sz: shape)
res*=sz;
return res;
}
cmav(const array<size_t,ndim> &shp_)
: shp(shp_), ptr(make_unique<vector<T>>(size())), d(ptr->data())
public:
mav_info(const shape_t &shape_, const stride_t &stride_)
: shp(shape_), str(stride_), sz(prod(shp))
{ MR_assert(shp.size()==str.size(), "dimensions mismatch"); }
mav_info(const shape_t &shape_)
: shp(shape_), sz(prod(shp))
{
MR_assert(ndim>0, "at least 1D required");
str[ndim-1]=1;
for (size_t i=2; i<=ndim; ++i)
str[ndim-i] = str[ndim-i+1]*shp[ndim-i+1];
}
cmav(const cmav &other) = default;
cmav(cmav &&other) = default;
operator cfmav<T>() const
{
return cfmav<T>(d, ptr, {shp.begin(), shp.end()}, {str.begin(), str.end()});
}
const T &operator[](size_t i) const
{ return operator()(i); }
const T &operator()(size_t i) const
{
static_assert(ndim==1, "ndim must be 1");
return d[str[0]*i];
}
const T &operator()(size_t i, size_t j) const
{
static_assert(ndim==2, "ndim must be 2");
return d[str[0]*i + str[1]*j];
}
const T &operator()(size_t i, size_t j, size_t k) const
{
static_assert(ndim==3, "ndim must be 3");
return d[str[0]*i + str[1]*j + str[2]*k];
str[ndim-i] = str[ndim-i+1]*ptrdiff_t(shp[ndim-i+1]);
}
// size_t ndim() const { return shp.size(); }
size_t size() const { return sz; }
const shape_t &shape() const { return shp; }
size_t shape(size_t i) const { return shp[i]; }
const array<size_t,ndim> &shape() const { return shp; }
size_t size() const
{
size_t res=1;
for (auto v: shp) res*=v;
return res;
}
ptrdiff_t stride(size_t i) const { return str[i]; }
const T *data() const
{ return d; }
const stride_t &stride() const { return str; }
const ptrdiff_t &stride(size_t i) const { return str[i]; }
bool last_contiguous() const
{ return (str[ndim-1]==1) || (str[ndim-1]==0); }
{ return (str.back()==1); }
bool contiguous() const
{
ptrdiff_t stride=1;
for (size_t i=0; i<ndim; ++i)
{
if (str[ndim-1-i]!=stride) return false;
stride *= shp[ndim-1-i];
stride *= ptrdiff_t(shp[ndim-1-i]);
}
return true;
}
template<typename T2> bool conformable(const cmav<T2,ndim> &other) const
bool conformable(const mav_info &other) const
{ return shp==other.shp; }
};
template<typename T, size_t ndim> class mav: public cmav<T, ndim>
template<typename T, size_t ndim> class mav: public mav_info<ndim>, public membuf<T>
{
static_assert((ndim>0) && (ndim<4), "only supports 1D, 2D, and 3D arrays");
using typename mav_info<ndim>::shape_t;
using typename mav_info<ndim>::stride_t;
protected:
using Tsp = shared_ptr<vector<T>>;
using parent = cmav<T, ndim>;
using parent::d;
using parent::ptr;
using parent::shp;
using parent::str;
using membuf<T>::Tsp;
using mav_info<ndim>::size;
using membuf<T>::d;
using membuf<T>::ptr;
using mav_info<ndim>::shp;
using mav_info<ndim>::str;
using membuf<T>::rw;
using membuf<T>::val;
public:
mav(T *d_, const array<size_t,ndim> &shp_,
const array<ptrdiff_t,ndim> &str_)
: parent(d_, shp_, str_) {}
mav(T *d_, const array<size_t,ndim> &shp_)
: parent(d_, shp_) {}
mav(const T *d_, const shape_t &shp_, const stride_t &str_)
: mav_info<ndim>(shp_, str_), membuf<T>(d_) {}
mav(T *d_, const shape_t &shp_, const stride_t &str_, bool rw_=false)
: mav_info<ndim>(shp_, str_), membuf<T>(d_, rw_) {}
mav(const T *d_, const shape_t &shp_)
: mav_info<ndim>(shp_), membuf<T>(d_) {}
mav(T *d_, const shape_t &shp_, bool rw_=false)
: mav_info<ndim>(shp_), membuf<T>(d_, rw_) {}
mav(const array<size_t,ndim> &shp_)
: parent(shp_) {}
mav(const mav &other) = default;
: mav_info<ndim>(shp_), membuf<T>(size()) {}
mav(const mav &other) = delete;
mav(mav &&other) = default;
operator fmav<T>() const
{
return fmav<T>(d, ptr, {shp.begin(), shp.end()}, {str.begin(), str.end()});
return fmav<T>(const_cast<T *>(d), ptr, {shp.begin(), shp.end()}, {str.begin(), str.end()}, rw);
}
T &operator[](size_t i) const
{ return operator()(i); }
T &operator()(size_t i) const
const T &operator()(size_t i) const
{
static_assert(ndim==1, "ndim must be 1");
return d[str[0]*i];
return val(str[0]*i);
}
T &operator()(size_t i, size_t j) const
const T &operator()(size_t i, size_t j) const