Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/core/static_types.hh
c74466647365features/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"
#include <type_traits>
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace detail {
template <UInt acc, UInt n, UInt... ns>
struct product_tail_rec : product_tail_rec<acc * n, ns...> {};
template <UInt acc, UInt n>
struct product_tail_rec<acc, n> : std::integral_constant<UInt, acc * n> {};
template <UInt N, UInt n, UInt... ns>
struct get_rec : get_rec<N - 1, ns...> {};
template <UInt n, UInt... ns>
struct get_rec<0, n, ns...> : std::integral_constant<UInt, n> {};
}
template <UInt... ns>
struct product : detail::product_tail_rec<1, ns...> {};
template <UInt N, UInt... ns>
struct get : detail::get_rec<N, ns...> {};
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;
using T_bare = typename std::remove_const<T>::type;
public:
using value_type = T;
static constexpr UInt size = _size;
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];
}
/// Scalar product
template <typename DT, typename ST>
__device__ __host__ T_bare dot(const StaticArray<DT, ST, size>& o) const {
decltype(T_bare(0) * DT(0)) res = 0;
for (UInt i = 0; i < size; ++i)
res += (*this)(i)*o(i);
return res;
}
/// L2 norm squared
__device__ __host__ T_bare l2squared() const { return this->dot(*this); }
/// L2 norm
__device__ __host__ T_bare 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 implicit copy operator
__device__ __host__ StaticArray& operator=(const StaticArray& o) {
return this->copy(o);
}
template <typename DT, typename ST>
__device__ __host__ void operator=(const StaticArray<DT, ST, size>& o) {
this->copy(o);
}
template <typename DT, typename ST>
__device__ __host__ StaticArray& copy(const StaticArray<DT, ST, size>& o) {
for (UInt i = 0; i < size; ++i)
(*this)(i) = o(i);
return *this;
}
T* begin() { return _mem; }
const T* begin() const { return _mem; }
T* end() { return _mem + size; }
const T* end() const { return _mem + size; }
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>
__device__ __host__ 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>
__device__ __host__ static UInt unpackOffset(UInt offset, UInt index) {
return offset + index;
}
public:
template <typename... Idx>
__device__ __host__ const T& operator()(Idx... idx) const {
return parent::operator()(unpackOffset(0, idx...));
}
template <typename... Idx>
__device__ __host__ T& operator()(Idx... idx) {
return parent::operator()(unpackOffset(0, idx...));
}
};
/* -------------------------------------------------------------------------- */
/* Common Static Types */
/* -------------------------------------------------------------------------- */
// Forward declaration
template <typename DataType, typename SupportType, UInt n>
class StaticVector;
template <typename DataType, typename SupportType, UInt n, UInt m>
class StaticMatrix : public StaticTensor<DataType, SupportType, n, m> {
using T = DataType;
using T_bare = typename std::remove_const<T>::type;
public:
using StaticTensor<DataType, SupportType, n, m>::operator=;
/// Outer product of two vectors
template <typename DT1, typename ST1, typename DT2, typename ST2>
__device__ __host__ void outer(const StaticVector<DT1, ST1, n>& a,
const StaticVector<DT2, ST2, m>& b);
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt l>
__device__ __host__ void mul(const StaticMatrix<DT1, ST1, n, l>& a,
const StaticMatrix<DT2, ST2, l, m>& b) {
(*this) = T(0);
for (UInt i = 0; i < n; ++i)
for (UInt j = 0; j < m; ++j)
for (UInt k = 0; k < l; ++k)
(*this)(i, j) += a(i, k) * b(k, j);
}
__device__ __host__ typename std::enable_if<n == m, T_bare>::type
trace() const {
T_bare res{0};
for (UInt i = 0 ; i < n; ++i)
res += (*this)(i, i);
return res;
}
template <typename DT1, typename ST1>
__device__ __host__ typename std::enable_if<n == m, void>::type
deviatoric(const StaticMatrix<DT1, ST1, n, m>& mat, Real factor = n) {
auto norm_trace = mat.trace() / factor;
for (UInt i = 0; i < n; ++i)
for (UInt j = 0; j < m; ++j)
(*this)(i, j) = mat(i, j) - (i == j) * norm_trace;
}
};
/// 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 <bool transpose, 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) +=
((transpose) ? mat(j, i) : mat(i, j)) * vec(j); // can be optimized
}
};
// Implementation of outer product
template <typename DataType, typename SupportType, UInt n, UInt m>
template <typename DT1, typename ST1, typename DT2, typename ST2>
__device__ __host__ void StaticMatrix<DataType, SupportType, n, m>::outer(
const StaticVector<DT1, ST1, n>& a, const StaticVector<DT2, ST2, m>& b) {
for (UInt i = 0; i < n; ++i)
for (UInt j = 0; j < m; ++j)
(*this)(i, j) = a(i) * b(j);
}
/* -------------------------------------------------------------------------- */
/* 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...> {
static constexpr UInt size = product<dims...>::value;
using parent = StaticParent<T, T[size], dims...>;
public:
using parent::operator=;
using parent::copy;
/// Default constructor
__device__ __host__ Tensor() = default;
/// Construct with default value
__device__ __host__ Tensor(T val) { *this = val; }
/// Construct from array
__device__ __host__ Tensor(const std::array<T, size>& arr) {
// we use size to ensure static loop unrolling
std::copy(arr.begin(), arr.begin() + size, this->_mem);
}
/// Copy from array
__device__ __host__ Tensor& operator=(const std::array<T, size>& arr) {
// we use size to ensure static loop unrolling
std::copy(arr.begin(), arr.begin() + size, this->_mem);
}
/// Construct by copy from static tensor
template <typename DT, typename ST>
__device__ __host__ 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>;
/* -------------------------------------------------------------------------- */
/* 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...> {
using parent = StaticParent<T, T*, dims...>;
public:
/// Implicit wrap constructor
__device__ __host__ explicit TensorProxy(T& spot) { this->_mem = &spot; }
/// Explicit construction from data location
__device__ __host__ explicit TensorProxy(T* spot) { this->_mem = spot; }
using parent::operator=;
public:
using stack_type = Tensor<StaticParent, T, dims...>;
};
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>;
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* Arithmetic operators creating temporaries */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* Simple operators */
/* -------------------------------------------------------------------------- */
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
__device__ __host__ Vector<decltype(DT1(0) + DT2(0)), dim>
operator+(const StaticVector<DT1, ST1, dim>& a,
const StaticVector<DT2, ST2, dim>& b) {
Vector<decltype(DT1(0) + DT2(0)), dim> res(a);
res += b;
return res;
}
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
__device__ __host__ Vector<decltype(DT1(0) - DT2(0)), dim>
operator-(const StaticVector<DT1, ST1, dim>& a,
const StaticVector<DT2, ST2, dim>& b) {
Vector<decltype(DT1(0) - DT2(0)), dim> res(a);
res -= b;
return res;
}
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
UInt m>
__device__ __host__ Matrix<decltype(DT1(0) + DT2(0)), n, m>
operator+(const StaticMatrix<DT1, ST1, n, m>& a,
const StaticMatrix<DT2, ST2, n, m>& b) {
Matrix<decltype(DT1(0) + DT2(0)), n, m> res(a);
res += b;
return res;
}
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
UInt m>
__device__ __host__ Matrix<decltype(DT1(0) - DT2(0)), n, m>
operator-(const StaticMatrix<DT1, ST1, n, m>& a,
const StaticMatrix<DT2, ST2, n, m>& b) {
Matrix<decltype(DT1(0) - DT2(0)), n, m> res(a);
res -= b;
return res;
}
template <
typename DT, typename ST, typename T, UInt n,
typename = typename std::enable_if<is_arithmetic<T>::value, void>::type>
Vector<decltype(DT(0) * T(0)), n> operator*(const StaticVector<DT, ST, n>& a,
const T& b) {
Vector<decltype(DT(0) * T(0)), n> res{a};
res *= b;
return res;
}
// symmetry
template <
typename DT, typename ST, typename T, UInt n,
typename = typename std::enable_if<is_arithmetic<T>::value, void>::type>
Vector<decltype(DT(0) * T(0)), n> operator*(const T& b,
const StaticVector<DT, ST, n>& a) {
return a * b;
}
template <
typename DT, typename ST, typename T, UInt n, UInt m,
typename = typename std::enable_if<is_arithmetic<T>::value, void>::type>
Matrix<decltype(DT(0) * T(0)), n, m>
operator*(const StaticMatrix<DT, ST, n, m>& a, const T& b) {
Matrix<decltype(DT(0) * T(0)), n, m> res{a};
res *= b;
return res;
}
// symmetry
template <
typename DT, typename ST, typename T, UInt n, UInt m,
typename = typename std::enable_if<is_arithmetic<T>::value, void>::type>
Matrix<decltype(DT(0) * T(0)), n, m>
operator*(const T& b, const StaticMatrix<DT, ST, n, m>& a) {
return a * b;
}
/* -------------------------------------------------------------------------- */
/* Linear algebra operators */
/* -------------------------------------------------------------------------- */
/// Matrix-vector multiplication
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
UInt m>
__device__ __host__ Vector<decltype(DT1(0) * DT2(0)), n>
operator*(const StaticMatrix<DT1, ST1, n, m>& a,
const StaticVector<DT2, ST2, m>& b) {
Vector<decltype(DT1(0) * DT2(0)), n> res;
res.template mul<false>(a, b);
return res;
}
/// Matrix-matrix multiplication
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
UInt m, UInt l>
__device__ __host__ Vector<decltype(DT1(0) * DT2(0)), n>
operator*(const StaticMatrix<DT1, ST1, n, l>& a,
const StaticMatrix<DT2, ST2, l, m>& b) {
Matrix<decltype(DT1(0) * DT2(0)), n, m> res;
res.template mul<false>(a, b);
return res;
}
template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
UInt m>
__device__ __host__ Matrix<decltype(DT1(0) * DT2(0)), n, m>
outer(const StaticVector<DT1, ST1, n>& a, const StaticVector<DT2, ST2, m>& b) {
Matrix<decltype(DT1(0) * DT2(0)), n, m> res;
res.outer(a, b);
return res;
}
__END_TAMAAS__
#endif // __STATIC_TYPES_HH__
c4science · Help