You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

404 lines
11 KiB

#ifndef ALGOIM_XARRAY_HPP
#define ALGOIM_XARRAY_HPP
// algoim::xarray implements an N-dimensional array view of a memory
// block controlled by the user
#include <cassert>
#include "uvector.hpp"
#include "sparkstack.hpp"
#include "utility.hpp"
namespace algoim
{
// MiniLoop<N> is an optimised version of MultiLoop<N> with zero lowerbound
template <int N>
class MiniLoop
{
uvector<int, N> i;
int iexp;
const uvector<int, N> ext;
public:
explicit MiniLoop(const uvector<int, N>& ext) : ext(ext), i(0), iexp(0) {}
MiniLoop& operator++()
{
++iexp;
for (int dim = N - 1; dim >= 0; --dim) {
if (++i(dim) < ext(dim)) return *this;
if (dim == 0) return *this;
i(dim) = 0;
}
return *this;
}
const uvector<int, N>& operator()() const { return i; }
int operator()(int index) const { return i(index); }
bool operator~() const { return i(0) < ext(0); }
int furl() const { return iexp; }
uvector<int, N> shifted(int dim, int amount) const
{
uvector<int, N> j = i;
j(dim) += amount;
return j;
}
};
template <typename T>
struct xarraySlice {
T* ptr;
int len;
xarraySlice(xarraySlice&) = delete;
xarraySlice(xarraySlice&&) = delete;
xarraySlice& operator=(const T& x)
{
for (int i = 0; i < len; ++i) ptr[i] = x;
return *this;
}
xarraySlice& operator+=(const T& x)
{
for (int i = 0; i < len; ++i) ptr[i] += x;
return *this;
}
xarraySlice& operator-=(const T& x)
{
for (int i = 0; i < len; ++i) ptr[i] -= x;
return *this;
}
xarraySlice& operator*=(const T& x)
{
for (int i = 0; i < len; ++i) ptr[i] *= x;
return *this;
}
xarraySlice& operator/=(const T& x)
{
for (int i = 0; i < len; ++i) ptr[i] /= x;
return *this;
}
xarraySlice& operator=(const xarraySlice& x)
{
for (int i = 0; i < len; ++i) ptr[i] = x.ptr[i];
return *this;
}
xarraySlice& operator+=(const xarraySlice& x)
{
for (int i = 0; i < len; ++i) ptr[i] += x.ptr[i];
return *this;
}
xarraySlice& operator-=(const xarraySlice& x)
{
for (int i = 0; i < len; ++i) ptr[i] -= x.ptr[i];
return *this;
}
xarraySlice& operator*=(const xarraySlice& x)
{
for (int i = 0; i < len; ++i) ptr[i] *= x.ptr[i];
return *this;
}
xarraySlice& operator/=(const xarraySlice& x)
{
for (int i = 0; i < len; ++i) ptr[i] /= x.ptr[i];
return *this;
}
template <typename X, typename Y>
struct prod {
const X& x;
const Y& y;
};
xarraySlice& operator=(const prod<xarraySlice, T>& op)
{
for (int i = 0; i < len; ++i) ptr[i] = op.x.ptr[i] * op.y;
return *this;
}
xarraySlice& operator+=(const prod<xarraySlice, T>& op)
{
for (int i = 0; i < len; ++i) ptr[i] += op.x.ptr[i] * op.y;
return *this;
}
xarraySlice& operator-=(const prod<xarraySlice, T>& op)
{
for (int i = 0; i < len; ++i) ptr[i] -= op.x.ptr[i] * op.y;
return *this;
}
xarraySlice& operator*=(const prod<xarraySlice, T>& op)
{
for (int i = 0; i < len; ++i) ptr[i] *= op.x.ptr[i] * op.y;
return *this;
}
xarraySlice& operator/=(const prod<xarraySlice, T>& op)
{
for (int i = 0; i < len; ++i) ptr[i] /= op.x.ptr[i] * op.y;
return *this;
}
};
template <typename T>
auto operator*(const xarraySlice<T>& x, const T& y)
{
return typename xarraySlice<T>::template prod<xarraySlice<T>, T>{x, y};
};
template <typename T>
void swap(const xarraySlice<T>& x, const xarraySlice<T>& y)
{
using std::swap;
for (int i = 0; i < x.len; ++i) swap(x.ptr[i], y.ptr[i]);
}
// algoim::xarray implements an N-dimensional array view of a memory
// block controlled by the user
template <typename T, int N>
class xarray
{
public:
T* data_;
uvector<int, N> ext_;
friend class SparkStack<T>;
// interpret the given block of memory as an N-dimensional array of the given extent
xarray(T* data, const uvector<int, N>& ext) : data_(data), ext_(ext) {}
xarray(const xarray& x)
{
// this->data_ = nullptr;
// this->ext_ = x.ext();
// if (x.data_) {
// algoim_spark_alloc(T, *this);
// *this = x;
// };
// 这里只能浅拷贝,因为不能在这里用algoim_spark_alloc分配新的内存
this->data_ = x.data_;
this->ext_ = x.ext_;
};
xarray()
{
this->data_ = nullptr;
this->ext_ = uvector<int, N>{};
}
xarray& operator=(const xarray& x)
{
// if (same_shape(x)) {
// if (x.data_) {
// for (int i = 0; i < size(); ++i) data_[i] = x.data_[i];
// }
// } else {
// this->ext_ = x.ext();
// if (x.data_) {
// algoim_spark_alloc(T, *this);
// *this = x;
// };
// }
assert(same_shape(x));
for (int i = 0; i < size(); ++i) data_[i] = x.data_[i];
return *this;
}
template <typename S>
xarray& operator=(const S& x)
{
for (int i = 0; i < size(); ++i) data_[i] = x;
return *this;
}
template <typename S>
xarray& operator+=(const S& x)
{
for (int i = 0; i < size(); ++i) data_[i] += x;
return *this;
}
xarray& operator+=(const xarray& x)
{
assert(same_shape(x));
for (int i = 0; i < size(); ++i) data_[i] += x.data_[i];
return *this;
}
template <typename S>
xarray& operator-=(const S& x)
{
for (int i = 0; i < size(); ++i) data_[i] -= x;
return *this;
}
xarray& operator-=(const xarray& x)
{
assert(same_shape(x));
for (int i = 0; i < size(); ++i) data_[i] -= x.data_[i];
return *this;
}
template <typename S>
xarray& operator*=(const S& x)
{
for (int i = 0; i < size(); ++i) data_[i] *= x;
return *this;
}
const T* data() const { return data_; }
T* data() { return data_; }
// Accessors for user-expanded index, 0 <= i < prod(ext)
const T& operator[](int i) const { return data_[i]; }
T& operator[](int i) { return data_[i]; }
// Slice and flatten this(i,:)
template <int NN = N, std::enable_if_t<NN == 1, bool> = true>
T& a(int i)
{
return data_[i];
}
// Slice and flatten this(i,:)
template <int NN = N, std::enable_if_t<NN == 1, bool> = true>
const T& a(int i) const
{
return data_[i];
}
// Slice and flatten this(i,:)
template <int NN = N, std::enable_if_t<(NN > 1), bool> = true>
auto a(int i) const
{
int span = prod(ext_, 0);
return xarraySlice<T>{data_ + i * span, span};
}
const T& operator()(int i, int j) const
{
static_assert(N == 2, "N == 2 required for integer pair access");
return data_[i * ext(1) + j];
}
T& operator()(int i, int j)
{
static_assert(N == 2, "N == 2 required for integer pair access");
return data_[i * ext(1) + j];
}
// Accessors by multi-index
const T& m(const uvector<int, N>& i) const { return data_[util::furl(i, ext_)]; }
// Accessors by multi-index
T& m(const uvector<int, N>& i) { return data_[util::furl(i, ext_)]; }
// Accessors by mini-loop, where it is assumed the MiniLoop object's
// associated extent is identical to this xarray
const T& l(const MiniLoop<N>& i) const { return data_[i.furl()]; }
T& l(const MiniLoop<N>& i) { return data_[i.furl()]; }
const uvector<int, N>& ext() const { return ext_; }
int ext(int i) const { return ext_(i); }
int size() const { return prod(ext_); }
MiniLoop<N> loop() const { return MiniLoop<N>(ext_); }
// xarray is strided so that the inner-most dimension (dim 0) is the slowest varying;
// as such, one may view an N-dimensionl xarray as a length ext(0) vector of (N-1)-
// dimensional xarrays. Two views of this(i,:) are provided:
// - flatten() returns a 2D view, of dimensions ext(0) by ext(1) * ... * ext(N - 1)
// - slice(i) returns an (N-1)-D array of extent (ext(1), ..., ext(N-1))
auto flatten() const { return xarray<T, 2>(data_, uvector<int, 2>{ext_(0), prod(ext_, 0)}); }
// xarray is strided so that the inner-most dimension (dim 0) is the slowest varying;
// as such, one may view an N-dimensionl xarray as a length ext(0) vector of (N-1)-
// dimensional xarrays. Two views of this(i,:) are provided:
// - flatten() returns a 2D view, of dimensions ext(0) by ext(1) * ... * ext(N - 1)
// - slice(i) returns an (N-1)-D array of extent (ext(1), ..., ext(N-1))
auto slice(int i) const { return xarray<T, N - 1>(data_ + i * prod(ext_, 0), remove_component(ext_, 0)); }
bool same_shape(const xarray& x) const { return all(x.ext_ == ext_); }
T maxNorm() const
{
assert(size() > 0);
using std::abs;
using std::max;
T m = abs(data_[0]);
for (int i = 1; i < size(); ++i) m = max(m, abs(data_[i]));
return m;
}
T min() const
{
assert(size() > 0);
using std::min;
T m = data_[0];
for (int i = 1; i < size(); ++i) m = min(m, data_[i]);
return m;
}
T max() const
{
assert(size() > 0);
using std::max;
T m = data_[0];
for (int i = 1; i < size(); ++i) m = max(m, data_[i]);
return m;
}
// Alter the extent of the view of the existing memory block; usually this
// only makes sense when the new extent is smaller than the existing extent,
// i.e., prod(new extent) <= prod(old extent)
void alterExtent(const uvector<int, N>& ext) { ext_ = ext; }
// flatten() and slice() return temporary objects which cannot be bound to non-const references;
// applying .ref() to the temporary object allows this to be done and is an assertion by the
// user that doing so is safe (relating to the output of flatten() or slice() being destroyed
// at the end of the full-expression)
xarray<T, N>& ref() { return *this; }
};
using tensor3 = xarray<real, 3>;
using tensor2 = xarray<real, 2>;
template <typename T, int N>
void xarrayInit(xarray<T, N>& arr)
{
for (auto i = arr.loop(); ~i; ++i) { arr.l(i) = 0; }
}
auto xarray2StdVector(const xarray<real, 3>& array)
{
auto ext = array.ext();
std::vector<std::vector<std::vector<real>>> v(ext(0), std::vector<std::vector<real>>(ext(1), std::vector<real>(ext(2), 0)));
for (int i = 0; i < ext(0); ++i) {
for (int j = 0; j < ext(1); ++j) {
for (int k = 0; k < ext(2); ++k) { v[i][j][k] = array.m(uvector<int, 3>(i, j, k)); }
}
}
return v;
}
} // namespace algoim
#endif