Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/core/static_types.hh
93e520e480b8features/netcdf-mpi

Last Change
Raw File

static_types.hh
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __STATIC_TYPES_HH__
#define __STATIC_TYPES_HH__
/* -------------------------------------------------------------------------- */
#include "tamaas.hh"
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace detail {
template <UInt acc, UInt n, UInt... ns>
struct product_tail_rec {
static constexpr UInt value = product_tail_rec<acc * n, ns...>::value;
};
template <UInt acc, UInt n>
struct product_tail_rec<acc, n> {
static constexpr UInt value = acc * n;
};
template <UInt N, UInt n, UInt... ns>
struct get_rec {
static constexpr UInt value = get_rec<N-1, ns...>::value;
};
template <UInt n, UInt... ns>
struct get_rec<0, n, ns...> {
static constexpr UInt value = n;
};
}
template <UInt... ns>
struct product {
static constexpr UInt value = detail::product_tail_rec<1, ns...>::value;
};
template <UInt N, UInt... ns>
struct get {
static constexpr UInt value = detail::get_rec<N, ns...>::value;
};
template <typename T>
struct is_arithmetic : std::is_arithmetic<T> {};
template <typename T>
struct is_arithmetic<thrust::complex<T>> : std::true_type {};
/**
* @brief Static Array
*
* This class is meant to be a small and fast object for intermediate
* calculations, possibly on wrapped memory belonging to a grid. Support type
* show be either a pointer or a C array. It should not contain any virtual
* method.
*/
template <typename DataType, typename SupportType, UInt size>
class StaticArray {
static_assert(std::is_array<SupportType>::value ||
std::is_pointer<SupportType>::value,
"the support type of StaticArray should be either a pointer or "
"a C-array");
using T = DataType;
public:
/// Access operator
__device__ __host__ T& operator()(UInt i) {
// TAMAAS_ASSERT(i < n, "Access out of bounds");
return _mem[i];
}
/// Access operator
__device__ __host__ const T& operator()(UInt i) const {
// TAMAAS_ASSERT(i < n, "Access out of bounds");
return _mem[i];
}
/// L2 norm squared
__device__ __host__ T l2squared() const {
T res = 0;
for (UInt i = 0; i < size; ++i)
res += _mem[i] * _mem[i];
return res;
}
/// L2 norm
__device__ __host__ T l2norm() const { return std::sqrt(l2squared()); }
#define VECTOR_OP(op) \
template <typename DT, typename ST> \
__device__ __host__ void operator op(const StaticArray<DT, ST, size>& o) { \
for (UInt i = 0; i < size; ++i) \
(*this)(i) op o(i); \
}
VECTOR_OP(+=)
VECTOR_OP(-=)
VECTOR_OP(*=)
VECTOR_OP(/=)
#undef VECTOR_OP
#define SCALAR_OP(op) \
template <typename T1> \
__device__ __host__ \
typename std::enable_if<is_arithmetic<T1>::value, StaticArray&>::type \
operator op(const T1& x) { \
for (UInt i = 0; i < size; ++i) \
(*this)(i) op x; \
return *this; \
}
SCALAR_OP(+=)
SCALAR_OP(-=)
SCALAR_OP(*=)
SCALAR_OP(/=)
SCALAR_OP(=)
#undef SCALAR_OP
/// Overriding the implit copy operator
StaticArray& operator=(const StaticArray& o) {
return this->copy(o);
}
template <typename DT, typename ST>
StaticArray& copy(const StaticArray<DT, ST, size>& o) {
for (UInt i = 0; i < size; ++i)
(*this)(i) = o(i);
return *this;
}
protected:
SupportType _mem;
};
/**
* @brief Static Tensor
*
* This class implements a multi-dimensional tensor behavior.
*/
template <typename DataType, typename SupportType = DataType*, UInt... dims>
class StaticTensor
: public StaticArray<DataType, SupportType, product<dims...>::value> {
using parent = StaticArray<DataType, SupportType, product<dims...>::value>;
using T = DataType;
public:
static constexpr UInt dim = sizeof...(dims);
using parent::operator=;
private:
template <typename... Idx>
static UInt unpackOffset(UInt offset, UInt index, Idx... rest) {
constexpr UInt size = sizeof...(rest);
offset += index;
offset *= get<dim-size, dims...>::value;
return unpackOffset(offset, rest...);
}
template <typename... Idx>
static UInt unpackOffset(UInt offset, UInt index) {
return offset + index;
}
public:
template <typename... Idx>
const T& operator()(Idx... idx) const {
return this->_mem[unpackOffset(0, idx...)];
}
template <typename... Idx>
T& operator()(Idx... idx) {
return this->_mem[unpackOffset(0, idx...)];
}
};
/* -------------------------------------------------------------------------- */
/* Common Static Types */
/* -------------------------------------------------------------------------- */
template <typename DataType, typename SupportType, UInt n, UInt m>
using StaticMatrix = StaticTensor<DataType, SupportType, n, m>;
/// Vector class with size determined at compile-time
template <typename DataType, typename SupportType, UInt n>
class StaticVector : public StaticTensor<DataType, SupportType, n> {
using T = DataType;
public:
using StaticTensor<DataType, SupportType, n>::operator=;
/// Matrix-vector product
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt m>
__device__ __host__ void mul(const StaticMatrix<DT1, ST1, n, m>& mat,
const StaticVector<DT2, ST2, m>& vec) {
*this = T(0);
for (UInt i = 0; i < n; ++i)
for (UInt j = 0; j < m; ++j)
(*this)(i) += mat(i, j) * vec(j);
}
};
/* -------------------------------------------------------------------------- */
/* Proxy Static Types */
/* -------------------------------------------------------------------------- */
/// Proxy type for tensor
template <template <typename, typename, UInt...> class StaticParent, typename T,
UInt... dims>
class TensorProxy : public StaticParent<T, T*, dims...> {
public:
/// Implicit wrap constructor
__device__ __host__ TensorProxy(T& spot) {
this->_mem = &spot;
}
using StaticTensor<T, T*, dims...>::operator=;
};
template <typename T, UInt n, UInt m>
using MatrixProxy = TensorProxy<StaticMatrix, T, n, m>;
template <typename T, UInt n>
using VectorProxy = TensorProxy<StaticVector, T, n>;
/* -------------------------------------------------------------------------- */
/* On the stack static types */
/* -------------------------------------------------------------------------- */
template <template <typename, typename, UInt...> class StaticParent, typename T,
UInt... dims>
class Tensor : public StaticParent<T, T[product<dims...>::value], dims...> {
using parent = StaticParent<T, T[product<dims...>::value], dims...>;
public:
using parent::operator=;
using parent::copy;
/// Default constructor
Tensor() = default;
/// Construct by copy from static tensor
template <typename DT, typename ST>
Tensor(const StaticParent<DT, ST, dims...> & o) {
this->copy(o);
}
};
template <typename T, UInt n, UInt m>
using Matrix = Tensor<StaticMatrix, T, n, m>;
template <typename T, UInt n>
using Vector = Tensor<StaticVector, T, n>;
/* -------------------------------------------------------------------------- */
/* Arithmetic operators creating temporaries */
/* -------------------------------------------------------------------------- */
template <template <typename, typename, UInt...> class StaticParent,
typename DT1, typename ST1, typename DT2, typename ST2, UInt... dims>
Tensor<StaticParent, decltype(DT1(0) + DT2(0)), dims...>
operator+(const StaticParent<DT1, ST1, dims...>& t1,
const StaticParent<DT2, ST2, dims...>& t2) {
Tensor<StaticParent, decltype(DT1(0) + DT2(0)), dims...> a(t1);
return a += t2;
}
__END_TAMAAS__
#endif // __STATIC_TYPES_HH__
c4science · Help