diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh index d194b5ca7..17d7f4207 100644 --- a/src/common/aka_types.hh +++ b/src/common/aka_types.hh @@ -1,1454 +1,1490 @@ /** * @file aka_types.hh * * @author Nicolas Richart * * @date creation: Thu Feb 17 2011 * @date last modification: Tue Feb 20 2018 * * @brief description of the "simple" types * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_error.hh" #include "aka_fwd.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_TYPES_HH__ #define __AKANTU_AKA_TYPES_HH__ namespace akantu { enum NormType { L_1 = 1, L_2 = 2, L_inf = UInt(-1) }; /** * DimHelper is a class to generalize the setup of a dim array from 3 * values. This gives a common interface in the TensorStorage class * independently of its derived inheritance (Vector, Matrix, Tensor3) * @tparam dim */ template struct DimHelper { static inline void setDims(UInt m, UInt n, UInt p, UInt dims[dim]); }; /* -------------------------------------------------------------------------- */ template <> struct DimHelper<1> { static inline void setDims(UInt m, __attribute__((unused)) UInt n, __attribute__((unused)) UInt p, UInt dims[1]) { dims[0] = m; } }; /* -------------------------------------------------------------------------- */ template <> struct DimHelper<2> { static inline void setDims(UInt m, UInt n, __attribute__((unused)) UInt p, UInt dims[2]) { dims[0] = m; dims[1] = n; } }; /* -------------------------------------------------------------------------- */ template <> struct DimHelper<3> { static inline void setDims(UInt m, UInt n, UInt p, UInt dims[3]) { dims[0] = m; dims[1] = n; dims[2] = p; } }; /* -------------------------------------------------------------------------- */ template class TensorStorage; /* -------------------------------------------------------------------------- */ /* Proxy classes */ /* -------------------------------------------------------------------------- */ namespace tensors { template struct is_copyable { enum : bool { value = false }; }; template struct is_copyable { enum : bool { value = true }; }; template struct is_copyable { enum : bool { value = true }; }; template struct is_copyable { enum : bool { value = true }; }; } // namespace tensors /** * @class TensorProxy aka_types.hh * @desc The TensorProxy class is a proxy class to the TensorStorage it handles * the * wrapped case. That is to say if an accessor should give access to a Tensor * wrapped on some data, like the Array::iterator they can return a * TensorProxy that will be automatically transformed as a TensorStorage wrapped * on the same data * @tparam T stored type * @tparam ndim order of the tensor * @tparam RetType real derived type */ template class TensorProxy { protected: using RetTypeProxy = typename _RetType::proxy; constexpr TensorProxy(T * data, UInt m, UInt n, UInt p) { DimHelper::setDims(m, n, p, this->n); this->values = data; } #ifndef SWIG template ::value>> explicit TensorProxy(const Other & other) { this->values = other.storage(); for (UInt i = 0; i < ndim; ++i) this->n[i] = other.size(i); } #endif public: using RetType = _RetType; UInt size(UInt i) const { AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim << " dimensions, not " << (i + 1)); return n[i]; } inline UInt size() const { UInt _size = 1; for (UInt d = 0; d < ndim; ++d) _size *= this->n[d]; return _size; } T * storage() const { return values; } #ifndef SWIG template ::value>> inline TensorProxy & operator=(const Other & other) { AKANTU_DEBUG_ASSERT( other.size() == this->size(), "You are trying to copy two tensors with different sizes"); memcpy(this->values, other.storage(), this->size() * sizeof(T)); return *this; } #endif // template ::value>> // inline TensorProxy & operator=(const Other && other) { // AKANTU_DEBUG_ASSERT( // other.size() == this->size(), // "You are trying to copy two tensors with different sizes"); // memcpy(this->values, other.storage(), this->size() * sizeof(T)); // return *this; // } template inline RetTypeProxy & operator*=(const O & o) { RetType(*this) *= o; return static_cast(*this); } template inline RetTypeProxy & operator/=(const O & o) { RetType(*this) /= o; return static_cast(*this); } protected: T * values; UInt n[ndim]; }; /* -------------------------------------------------------------------------- */ template class VectorProxy : public TensorProxy> { using parent = TensorProxy>; using type = Vector; public: constexpr VectorProxy(T * data, UInt n) : parent(data, n, 0, 0) {} template explicit VectorProxy(Other & src) : parent(src) {} /* ---------------------------------------------------------------------- */ template inline VectorProxy & operator=(const Other & other) { parent::operator=(other); return *this; } // inline VectorProxy & operator=(const VectorProxy && other) { // parent::operator=(other); // return *this; // } /* ------------------------------------------------------------------------ */ T & operator()(UInt index) { return this->values[index]; }; const T & operator()(UInt index) const { return this->values[index]; }; }; template class MatrixProxy : public TensorProxy> { using parent = TensorProxy>; using type = Matrix; public: MatrixProxy(T * data, UInt m, UInt n) : parent(data, m, n, 0) {} template explicit MatrixProxy(Other & src) : parent(src) {} /* ---------------------------------------------------------------------- */ template inline MatrixProxy & operator=(const Other & other) { parent::operator=(other); return *this; } }; template class Tensor3Proxy : public TensorProxy> { using parent = TensorProxy>; using type = Tensor3; public: Tensor3Proxy(const T * data, UInt m, UInt n, UInt k) : parent(data, m, n, k) {} Tensor3Proxy(const Tensor3Proxy & src) : parent(src) {} Tensor3Proxy(const Tensor3 & src) : parent(src) {} /* ---------------------------------------------------------------------- */ template inline Tensor3Proxy & operator=(const Other & other) { parent::operator=(other); return *this; } }; /* -------------------------------------------------------------------------- */ /* Tensor base class */ /* -------------------------------------------------------------------------- */ template class TensorStorage : public TensorTrait { public: using value_type = T; friend class Array; protected: template void copySize(const TensorType & src) { for (UInt d = 0; d < ndim; ++d) this->n[d] = src.size(d); this->_size = src.size(); } TensorStorage() : values(nullptr) { for (UInt d = 0; d < ndim; ++d) this->n[d] = 0; _size = 0; } TensorStorage(const TensorProxy & proxy) { this->copySize(proxy); this->values = proxy.storage(); this->wrapped = true; } public: TensorStorage(const TensorStorage & src) = delete; TensorStorage(const TensorStorage & src, bool deep_copy) : values(nullptr) { if (deep_copy) this->deepCopy(src); else this->shallowCopy(src); } protected: TensorStorage(UInt m, UInt n, UInt p, const T & def) { static_assert(std::is_trivially_constructible{}, "Cannot create a tensor on non trivial types"); DimHelper::setDims(m, n, p, this->n); this->computeSize(); this->values = new T[this->_size]; this->set(def); this->wrapped = false; } TensorStorage(T * data, UInt m, UInt n, UInt p) { DimHelper::setDims(m, n, p, this->n); this->computeSize(); this->values = data; this->wrapped = true; } public: /* ------------------------------------------------------------------------ */ template inline void shallowCopy(const TensorType & src) { this->copySize(src); if (!this->wrapped) delete[] this->values; this->values = src.storage(); this->wrapped = true; } /* ------------------------------------------------------------------------ */ template inline void deepCopy(const TensorType & src) { this->copySize(src); if (!this->wrapped) delete[] this->values; static_assert(std::is_trivially_constructible{}, "Cannot create a tensor on non trivial types"); this->values = new T[this->_size]; static_assert(std::is_trivially_copyable{}, "Cannot copy a tensor on non trivial types"); memcpy((void *)this->values, (void *)src.storage(), this->_size * sizeof(T)); this->wrapped = false; } virtual ~TensorStorage() { if (!this->wrapped) delete[] this->values; } /* ------------------------------------------------------------------------ */ inline TensorStorage & operator=(const TensorStorage & src) { return this->operator=(dynamic_cast(src)); } /* ------------------------------------------------------------------------ */ inline TensorStorage & operator=(const RetType & src) { if (this != &src) { if (this->wrapped) { static_assert(std::is_trivially_copyable{}, "Cannot copy a tensor on non trivial types"); // this test is not sufficient for Tensor of order higher than 1 AKANTU_DEBUG_ASSERT(this->_size == src.size(), "Tensors of different size"); memcpy((void *)this->values, (void *)src.storage(), this->_size * sizeof(T)); } else { deepCopy(src); } } return *this; } /* ------------------------------------------------------------------------ */ template inline RetType & operator+=(const TensorStorage & other) { T * a = this->storage(); T * b = other.storage(); AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); for (UInt i = 0; i < _size; ++i) *(a++) += *(b++); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ template inline RetType & operator-=(const TensorStorage & other) { T * a = this->storage(); T * b = other.storage(); AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); for (UInt i = 0; i < _size; ++i) *(a++) -= *(b++); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator+=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) += x; return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator-=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) -= x; return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator*=(const T & x) { T * a = this->storage(); for (UInt i = 0; i < _size; ++i) *(a++) *= x; return *(static_cast(this)); } /* ---------------------------------------------------------------------- */ inline RetType & operator/=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) /= x; return *(static_cast(this)); } /// Y = \alpha X + Y inline RetType & aXplusY(const TensorStorage & other, const T & alpha = 1.) { AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); Math::aXplusY(this->_size, alpha, other.storage(), this->storage()); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ T * storage() const { return values; } UInt size() const { return _size; } UInt size(UInt i) const { AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim << " dimensions, not " << (i + 1)); return n[i]; }; /* ------------------------------------------------------------------------ */ inline void clear() { memset(values, 0, _size * sizeof(T)); }; inline void set(const T & t) { std::fill_n(values, _size, t); }; template inline void copy(const TensorType & other) { AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be copied"); memcpy(values, other.storage(), _size * sizeof(T)); } bool isWrapped() const { return this->wrapped; } protected: inline void computeSize() { _size = 1; for (UInt d = 0; d < ndim; ++d) _size *= this->n[d]; } protected: template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += std::pow(std::abs(*it), norm_type); return std::pow(_norm, 1. / norm_type); } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += std::abs(*it); return _norm; } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += *it * *it; return sqrt(_norm); } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm = std::max(std::abs(*it), _norm); return _norm; } }; public: /*----------------------------------------------------------------------- */ /// "Entrywise" norm norm @f[ \|\boldsymbol{T}\|_p = \left( /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}} /// @f] template inline T norm() const { return NormHelper::norm(*this); } protected: UInt n[ndim]; UInt _size; T * values; bool wrapped{false}; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ namespace types { namespace details { template class vector_iterator { public: using difference_type = std::ptrdiff_t; using value_type = std::decay_t; using pointer = value_type *; using reference = reference_; using iterator_category = std::input_iterator_tag; vector_iterator(pointer ptr) : ptr(ptr) {} // input iterator ++it vector_iterator & operator++() { ++ptr; return *this; } // input iterator it++ vector_iterator operator++(int) { auto cpy = *this; ++ptr; return cpy; } vector_iterator & operator+=(int n) { ptr += n; return *this; } vector_iterator operator+(int n) { vector_iterator cpy(*this); cpy += n; return cpy; } // input iterator it != other_it bool operator!=(const vector_iterator & other) const { return ptr != other.ptr; } bool operator==(const vector_iterator & other) const { return ptr == other.ptr; } difference_type operator-(const vector_iterator & other) const { return this->ptr - other.ptr; } // input iterator dereference *it reference operator*() { return *ptr; } pointer operator->() { return ptr; } private: pointer ptr; }; } // namespace details } // namespace types /* -------------------------------------------------------------------------- */ /* Vector */ /* -------------------------------------------------------------------------- */ template class Vector : public TensorStorage> { using parent = TensorStorage>; public: using value_type = typename parent::value_type; using proxy = VectorProxy; public: Vector() : parent() {} explicit Vector(UInt n, const T & def = T()) : parent(n, 0, 0, def) {} Vector(T * data, UInt n) : parent(data, n, 0, 0) {} Vector(const Vector & src, bool deep_copy = true) : parent(src, deep_copy) {} Vector(const TensorProxy & src) : parent(src) {} Vector(std::initializer_list list) : parent(list.size(), 0, 0, T()) { UInt i = 0; for (auto val : list) { operator()(i++) = val; } } public: using iterator = types::details::vector_iterator; using const_iterator = types::details::vector_iterator; iterator begin() { return iterator(this->storage()); } iterator end() { return iterator(this->storage() + this->size()); } const_iterator begin() const { return const_iterator(this->storage()); } const_iterator end() const { return const_iterator(this->storage() + this->size()); } public: ~Vector() override = default; /* ------------------------------------------------------------------------ */ inline Vector & operator=(const Vector & src) { parent::operator=(src); return *this; } /* ------------------------------------------------------------------------ */ inline T & operator()(UInt i) { AKANTU_DEBUG_ASSERT((i < this->n[0]), "Access out of the vector! " << "Index (" << i << ") is out of the vector of size (" << this->n[0] << ")"); return *(this->values + i); } inline const T & operator()(UInt i) const { AKANTU_DEBUG_ASSERT((i < this->n[0]), "Access out of the vector! " << "Index (" << i << ") is out of the vector of size (" << this->n[0] << ")"); return *(this->values + i); } inline T & operator[](UInt i) { return this->operator()(i); } inline const T & operator[](UInt i) const { return this->operator()(i); } /* ------------------------------------------------------------------------ */ inline Vector & operator*=(Real x) { return parent::operator*=(x); } inline Vector & operator/=(Real x) { return parent::operator/=(x); } /* ------------------------------------------------------------------------ */ inline Vector & operator*=(const Vector & vect) { AKANTU_DEBUG_ASSERT(this->_size == vect._size, "The vectors have non matching sizes"); T * a = this->storage(); T * b = vect.storage(); for (UInt i = 0; i < this->_size; ++i) *(a++) *= *(b++); return *this; } /* ------------------------------------------------------------------------ */ inline Real dot(const Vector & vect) const { return Math::vectorDot(this->values, vect.storage(), this->_size); } /* ------------------------------------------------------------------------ */ inline Real mean() const { Real mean = 0; T * a = this->storage(); for (UInt i = 0; i < this->_size; ++i) mean += *(a++); return mean / this->_size; } /* ------------------------------------------------------------------------ */ inline Vector & crossProduct(const Vector & v1, const Vector & v2) { AKANTU_DEBUG_ASSERT(this->size() == 3, "crossProduct is only defined in 3D (n=" << this->size() << ")"); AKANTU_DEBUG_ASSERT( this->size() == v1.size() && this->size() == v2.size(), "crossProduct is not a valid operation non matching size vectors"); Math::vectorProduct3(v1.storage(), v2.storage(), this->values); return *this; } inline Vector crossProduct(const Vector & v) { Vector tmp(this->size()); tmp.crossProduct(*this, v); return tmp; } /* ------------------------------------------------------------------------ */ inline void solve(const Matrix & A, const Vector & b) { AKANTU_DEBUG_ASSERT( this->size() == A.rows() && this->_size == A.cols(), "The size of the solution vector mismatches the size of the matrix"); AKANTU_DEBUG_ASSERT( this->_size == b._size, "The rhs vector has a mismatch in size with the matrix"); Math::solve(this->_size, A.storage(), this->values, b.storage()); } /* ------------------------------------------------------------------------ */ template inline void mul(const Matrix & A, const Vector & x, Real alpha = 1.0); /* ------------------------------------------------------------------------ */ inline Real norm() const { return parent::template norm(); } template inline Real norm() const { return parent::template norm(); } /* ------------------------------------------------------------------------ */ inline Vector & normalize() { Real n = norm(); operator/=(n); return *this; } /* ------------------------------------------------------------------------ */ /// norm of (*this - x) inline Real distance(const Vector & y) const { Real * vx = this->values; Real * vy = y.storage(); Real sum_2 = 0; for (UInt i = 0; i < this->_size; ++i, ++vx, ++vy) sum_2 += (*vx - *vy) * (*vx - *vy); return sqrt(sum_2); } /* ------------------------------------------------------------------------ */ inline bool equal(const Vector & v, Real tolerance = Math::getTolerance()) const { T * a = this->storage(); T * b = v.storage(); UInt i = 0; while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance)) ++i; return i == this->_size; } /* ------------------------------------------------------------------------ */ inline short compare(const Vector & v, Real tolerance = Math::getTolerance()) const { T * a = this->storage(); T * b = v.storage(); for (UInt i(0); i < this->_size; ++i, ++a, ++b) { if (std::abs(*a - *b) > tolerance) return (((*a - *b) > tolerance) ? 1 : -1); } return 0; } /* ------------------------------------------------------------------------ */ inline bool operator==(const Vector & v) const { return equal(v); } inline bool operator!=(const Vector & v) const { return !operator==(v); } inline bool operator<(const Vector & v) const { return compare(v) == -1; } inline bool operator>(const Vector & v) const { return compare(v) == 1; } #ifndef SWIG template decltype(auto) accumulate(const Vector & v, Acc && accumulator, Func && func) const { T * a = this->storage(); T * b = v.storage(); for (UInt i(0); i < this->_size; ++i, ++a, ++b) { accumulator = func(*a, *b, std::forward(accumulator)); } return accumulator; } inline bool operator<=(const Vector & v) const { bool res = true; return accumulate(v, res, [](auto && a, auto && b, auto && accumulator) { return accumulator & (a <= b); }); } inline bool operator>=(const Vector & v) const { bool res = true; return accumulate(v, res, [](auto && a, auto && b, auto && accumulator) { return accumulator & (a >= b); }); } #endif /* ------------------------------------------------------------------------ */ /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << "["; for (UInt i = 0; i < this->_size; ++i) { if (i != 0) stream << ", "; stream << this->values[i]; } stream << "]"; } /* ---------------------------------------------------------------------- */ static inline Vector zeros(UInt n) { Vector tmp(n); tmp.set(T()); return tmp; } }; using RVector = Vector; /* ------------------------------------------------------------------------ */ template <> inline bool Vector::equal(const Vector & v, __attribute__((unused)) Real tolerance) const { UInt * a = this->storage(); UInt * b = v.storage(); UInt i = 0; while (i < this->_size && (*(a++) == *(b++))) ++i; return i == this->_size; } +/* -------------------------------------------------------------------------- */ +namespace types { + namespace details { + template class column_iterator { + public: + using difference_type = std::ptrdiff_t; + using value_type = decltype(std::declval().operator()(0)); + using pointer = value_type *; + using reference = value_type &; + using iterator_category = std::input_iterator_tag; + + + column_iterator(Mat & mat, UInt col) : mat(mat), col(col) {} + decltype(auto) operator*() { return mat(col); } + decltype(auto) operator++() { + ++col; + AKANTU_DEBUG_ASSERT(col <= mat.cols(), "The iterator is out of bound"); + return *this; + } + decltype(auto) operator++(int) { + auto tmp = *this; + ++col; + AKANTU_DEBUG_ASSERT(col <= mat.cols(), "The iterator is out of bound"); + return tmp; + } + bool operator!=(const column_iterator & other) const { + return col != other.col; + } + + bool operator==(const column_iterator & other) const { + return not operator!=(other); + } + + private: + Mat & mat; + UInt col; + }; + } // namespace details +} // namespace types + /* ------------------------------------------------------------------------ */ /* Matrix */ /* ------------------------------------------------------------------------ */ template class Matrix : public TensorStorage> { using parent = TensorStorage>; public: using value_type = typename parent::value_type; using proxy = MatrixProxy; public: Matrix() : parent() {} Matrix(UInt m, UInt n, const T & def = T()) : parent(m, n, 0, def) {} Matrix(T * data, UInt m, UInt n) : parent(data, m, n, 0) {} Matrix(const Matrix & src, bool deep_copy = true) : parent(src, deep_copy) {} Matrix(const MatrixProxy & src) : parent(src) {} Matrix(std::initializer_list> list) { static_assert(std::is_trivially_copyable{}, "Cannot create a tensor on non trivial types"); std::size_t n = 0; std::size_t m = list.size(); for (auto row : list) { n = std::max(n, row.size()); } DimHelper<2>::setDims(m, n, 0, this->n); this->computeSize(); this->values = new T[this->_size]; this->set(0); UInt i = 0, j = 0; for (auto & row : list) { for (auto & val : row) { at(i, j++) = val; } ++i; j = 0; } } ~Matrix() override = default; /* ------------------------------------------------------------------------ */ inline Matrix & operator=(const Matrix & src) { parent::operator=(src); return *this; } public: /* ---------------------------------------------------------------------- */ UInt rows() const { return this->n[0]; } UInt cols() const { return this->n[1]; } /* ---------------------------------------------------------------------- */ inline T & at(UInt i, UInt j) { AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])), "Access out of the matrix! " << "Index (" << i << ", " << j << ") is out of the matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return *(this->values + i + j * this->n[0]); } inline const T & at(UInt i, UInt j) const { AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])), "Access out of the matrix! " << "Index (" << i << ", " << j << ") is out of the matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return *(this->values + i + j * this->n[0]); } /* ------------------------------------------------------------------------ */ inline T & operator()(UInt i, UInt j) { return this->at(i, j); } inline const T & operator()(UInt i, UInt j) const { return this->at(i, j); } /// give a line vector wrapped on the column i inline VectorProxy operator()(UInt j) { AKANTU_DEBUG_ASSERT(j < this->n[1], "Access out of the matrix! " << "You are trying to access the column vector " << j << " in a matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return VectorProxy(this->values + j * this->n[0], this->n[0]); } inline const VectorProxy operator()(UInt j) const { AKANTU_DEBUG_ASSERT(j < this->n[1], "Access out of the matrix! " << "You are trying to access the column vector " << j << " in a matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return VectorProxy(this->values + j * this->n[0], this->n[0]); } #ifndef SWIG -private: - template class col_iterator { - public: - col_iterator(Mat & mat, UInt col) : mat(mat), col(col) {} - decltype(auto) operator*() { return mat(col); } - decltype(auto) operator++() { - ++col; - AKANTU_DEBUG_ASSERT(col <= mat.cols(), "The iterator is out of bound"); - return *this; - } - decltype(auto) operator++(int) { - auto tmp = *this; - ++col; - AKANTU_DEBUG_ASSERT(col <= mat.cols(), "The iterator is out of bound"); - return tmp; - } - bool operator!=(const col_iterator & other) { return col != other.col; } - - private: - Mat & mat; - UInt col; - }; - public: - decltype(auto) begin() { return col_iterator>(*this, 0); } + decltype(auto) begin() { + return types::details::column_iterator>(*this, 0); + } decltype(auto) begin() const { - return col_iterator>(*this, 0); + return types::details::column_iterator>(*this, 0); } - decltype(auto) end() { return col_iterator>(*this, this->cols()); } + decltype(auto) end() { + return types::details::column_iterator>(*this, this->cols()); + } decltype(auto) end() const { - return col_iterator>(*this, this->cols()); + return types::details::column_iterator>(*this, + this->cols()); } #endif /* ------------------------------------------------------------------------ */ inline void block(const Matrix & block, UInt pos_i, UInt pos_j) { AKANTU_DEBUG_ASSERT(pos_i + block.rows() <= rows(), "The block size or position are not correct"); AKANTU_DEBUG_ASSERT(pos_i + block.cols() <= cols(), "The block size or position are not correct"); for (UInt i = 0; i < block.rows(); ++i) for (UInt j = 0; j < block.cols(); ++j) this->at(i + pos_i, j + pos_j) = block(i, j); } inline Matrix block(UInt pos_i, UInt pos_j, UInt block_rows, UInt block_cols) const { AKANTU_DEBUG_ASSERT(pos_i + block_rows <= rows(), "The block size or position are not correct"); AKANTU_DEBUG_ASSERT(pos_i + block_cols <= cols(), "The block size or position are not correct"); Matrix block(block_rows, block_cols); for (UInt i = 0; i < block_rows; ++i) for (UInt j = 0; j < block_cols; ++j) block(i, j) = this->at(i + pos_i, j + pos_j); return block; } inline T & operator[](UInt idx) { return *(this->values + idx); }; inline const T & operator[](UInt idx) const { return *(this->values + idx); }; /* ---------------------------------------------------------------------- */ inline Matrix operator*(const Matrix & B) { Matrix C(this->rows(), B.cols()); C.mul(*this, B); return C; } /* ----------------------------------------------------------------------- */ inline Matrix & operator*=(const T & x) { return parent::operator*=(x); } inline Matrix & operator*=(const Matrix & B) { Matrix C(*this); this->mul(C, B); return *this; } /* ---------------------------------------------------------------------- */ template inline void mul(const Matrix & A, const Matrix & B, T alpha = 1.0) { UInt k = A.cols(); if (tr_A) k = A.rows(); #ifndef AKANTU_NDEBUG if (tr_B) { AKANTU_DEBUG_ASSERT(k == B.cols(), "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->cols() == B.rows(), "matrices to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(k == B.rows(), "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->cols() == B.cols(), "matrices to multiply have no fit dimensions"); } if (tr_A) { AKANTU_DEBUG_ASSERT(this->rows() == A.cols(), "matrices to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(this->rows() == A.rows(), "matrices to multiply have no fit dimensions"); } #endif // AKANTU_NDEBUG Math::matMul(this->rows(), this->cols(), k, alpha, A.storage(), B.storage(), 0., this->storage()); } /* ---------------------------------------------------------------------- */ inline void outerProduct(const Vector & A, const Vector & B) { AKANTU_DEBUG_ASSERT( A.size() == this->rows() && B.size() == this->cols(), "A and B are not compatible with the size of the matrix"); for (UInt i = 0; i < this->rows(); ++i) { for (UInt j = 0; j < this->cols(); ++j) { this->values[i + j * this->rows()] += A[i] * B[j]; } } } private: class EigenSorter { public: EigenSorter(const Vector & eigs) : eigs(eigs) {} bool operator()(const UInt & a, const UInt & b) const { return (eigs(a) > eigs(b)); } private: const Vector & eigs; }; public: /* ---------------------------------------------------------------------- */ inline void eig(Vector & eigenvalues, Matrix & eigenvectors) const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "eig is not a valid operation on a rectangular matrix"); AKANTU_DEBUG_ASSERT(eigenvalues.size() == this->cols(), "eigenvalues should be of size " << this->cols() << "."); #ifndef AKANTU_NDEBUG if (eigenvectors.storage() != nullptr) AKANTU_DEBUG_ASSERT((eigenvectors.rows() == eigenvectors.cols()) && (eigenvectors.rows() == this->cols()), "Eigenvectors needs to be a square matrix of size " << this->cols() << " x " << this->cols() << "."); #endif Matrix tmp = *this; Vector tmp_eigs(eigenvalues.size()); Matrix tmp_eig_vects(eigenvectors.rows(), eigenvectors.cols()); if (tmp_eig_vects.rows() == 0 || tmp_eig_vects.cols() == 0) Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage()); else Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage(), tmp_eig_vects.storage()); Vector perm(eigenvalues.size()); for (UInt i = 0; i < perm.size(); ++i) perm(i) = i; std::sort(perm.storage(), perm.storage() + perm.size(), EigenSorter(tmp_eigs)); for (UInt i = 0; i < perm.size(); ++i) eigenvalues(i) = tmp_eigs(perm(i)); if (tmp_eig_vects.rows() != 0 && tmp_eig_vects.cols() != 0) for (UInt i = 0; i < perm.size(); ++i) { for (UInt j = 0; j < eigenvectors.rows(); ++j) { eigenvectors(j, i) = tmp_eig_vects(j, perm(i)); } } } /* ---------------------------------------------------------------------- */ inline void eig(Vector & eigenvalues) const { Matrix empty; eig(eigenvalues, empty); } /* ---------------------------------------------------------------------- */ inline void eye(T alpha = 1.) { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "eye is not a valid operation on a rectangular matrix"); this->clear(); for (UInt i = 0; i < this->cols(); ++i) { this->values[i + i * this->rows()] = alpha; } } /* ---------------------------------------------------------------------- */ static inline Matrix eye(UInt m, T alpha = 1.) { Matrix tmp(m, m); tmp.eye(alpha); return tmp; } /* ---------------------------------------------------------------------- */ inline T trace() const { AKANTU_DEBUG_ASSERT( this->cols() == this->rows(), "trace is not a valid operation on a rectangular matrix"); T trace = 0.; for (UInt i = 0; i < this->rows(); ++i) { trace += this->values[i + i * this->rows()]; } return trace; } /* ---------------------------------------------------------------------- */ inline Matrix transpose() const { Matrix tmp(this->cols(), this->rows()); for (UInt i = 0; i < this->rows(); ++i) { for (UInt j = 0; j < this->cols(); ++j) { tmp(j, i) = operator()(i, j); } } return tmp; } /* ---------------------------------------------------------------------- */ inline void inverse(const Matrix & A) { AKANTU_DEBUG_ASSERT(A.cols() == A.rows(), "inv is not a valid operation on a rectangular matrix"); AKANTU_DEBUG_ASSERT(this->cols() == A.cols(), "the matrix should have the same size as its inverse"); if (this->cols() == 1) *this->values = 1. / *A.storage(); else if (this->cols() == 2) Math::inv2(A.storage(), this->values); else if (this->cols() == 3) Math::inv3(A.storage(), this->values); else Math::inv(this->cols(), A.storage(), this->values); } inline Matrix inverse() { Matrix inv(this->rows(), this->cols()); inv.inverse(*this); return inv; } /* --------------------------------------------------------------------- */ inline T det() const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "inv is not a valid operation on a rectangular matrix"); if (this->cols() == 1) return *(this->values); else if (this->cols() == 2) return Math::det2(this->values); else if (this->cols() == 3) return Math::det3(this->values); else return Math::det(this->cols(), this->values); } /* --------------------------------------------------------------------- */ inline T doubleDot(const Matrix & other) const { AKANTU_DEBUG_ASSERT( this->cols() == this->rows(), "doubleDot is not a valid operation on a rectangular matrix"); if (this->cols() == 1) return *(this->values) * *(other.storage()); else if (this->cols() == 2) return Math::matrixDoubleDot22(this->values, other.storage()); else if (this->cols() == 3) return Math::matrixDoubleDot33(this->values, other.storage()); else AKANTU_ERROR("doubleDot is not defined for other spatial dimensions" << " than 1, 2 or 3."); return T(); } /* ---------------------------------------------------------------------- */ /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << "["; for (UInt i = 0; i < this->n[0]; ++i) { if (i != 0) stream << ", "; stream << "["; for (UInt j = 0; j < this->n[1]; ++j) { if (j != 0) stream << ", "; stream << operator()(i, j); } stream << "]"; } stream << "]"; }; }; /* ------------------------------------------------------------------------ */ template template inline void Vector::mul(const Matrix & A, const Vector & x, Real alpha) { #ifndef AKANTU_NDEBUG UInt n = x.size(); if (tr_A) { AKANTU_DEBUG_ASSERT(n == A.rows(), "matrix and vector to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->size() == A.cols(), "matrix and vector to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(n == A.cols(), "matrix and vector to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->size() == A.rows(), "matrix and vector to multiply have no fit dimensions"); } #endif Math::matVectMul(A.rows(), A.cols(), alpha, A.storage(), x.storage(), 0., this->storage()); } /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Matrix & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Vector & _this) { _this.printself(stream); return stream; } /* ------------------------------------------------------------------------ */ /* Tensor3 */ /* ------------------------------------------------------------------------ */ template class Tensor3 : public TensorStorage> { using parent = TensorStorage>; public: using value_type = typename parent::value_type; using proxy = Tensor3Proxy; public: Tensor3() : parent(){}; Tensor3(UInt m, UInt n, UInt p, const T & def = T()) : parent(m, n, p, def) {} Tensor3(T * data, UInt m, UInt n, UInt p) : parent(data, m, n, p) {} Tensor3(const Tensor3 & src, bool deep_copy = true) : parent(src, deep_copy) {} Tensor3(const proxy & src) : parent(src) {} public: /* ------------------------------------------------------------------------ */ inline Tensor3 & operator=(const Tensor3 & src) { parent::operator=(src); return *this; } /* ---------------------------------------------------------------------- */ inline T & operator()(UInt i, UInt j, UInt k) { AKANTU_DEBUG_ASSERT( (i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the element " << "(" << i << ", " << j << ", " << k << ") in a tensor of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return *(this->values + (k * this->n[0] + i) * this->n[1] + j); } inline const T & operator()(UInt i, UInt j, UInt k) const { AKANTU_DEBUG_ASSERT( (i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the element " << "(" << i << ", " << j << ", " << k << ") in a tensor of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return *(this->values + (k * this->n[0] + i) * this->n[1] + j); } inline MatrixProxy operator()(UInt k) { AKANTU_DEBUG_ASSERT((k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the slice " << k << " in a tensor3 of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } inline const MatrixProxy operator()(UInt k) const { AKANTU_DEBUG_ASSERT((k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the slice " << k << " in a tensor3 of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } inline MatrixProxy operator[](UInt k) { return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } inline const MatrixProxy operator[](UInt k) const { return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } }; /* -------------------------------------------------------------------------- */ // support operations for the creation of other vectors /* -------------------------------------------------------------------------- */ template Vector operator*(const T & scalar, const Vector & a) { Vector r(a); r *= scalar; return r; } template Vector operator*(const Vector & a, const T & scalar) { Vector r(a); r *= scalar; return r; } template Vector operator/(const Vector & a, const T & scalar) { Vector r(a); r /= scalar; return r; } template Vector operator*(const Vector & a, const Vector & b) { Vector r(a); r *= b; return r; } template Vector operator+(const Vector & a, const Vector & b) { Vector r(a); r += b; return r; } template Vector operator-(const Vector & a, const Vector & b) { Vector r(a); r -= b; return r; } template Vector operator*(const Matrix & A, const Vector & b) { Vector r(b.size()); r.template mul(A, b); return r; } /* -------------------------------------------------------------------------- */ template Matrix operator*(const T & scalar, const Matrix & a) { Matrix r(a); r *= scalar; return r; } template Matrix operator*(const Matrix & a, const T & scalar) { Matrix r(a); r *= scalar; return r; } template Matrix operator/(const Matrix & a, const T & scalar) { Matrix r(a); r /= scalar; return r; } template Matrix operator+(const Matrix & a, const Matrix & b) { Matrix r(a); r += b; return r; } template Matrix operator-(const Matrix & a, const Matrix & b) { Matrix r(a); r -= b; return r; } } // namespace akantu #include namespace std { template struct iterator_traits<::akantu::types::details::vector_iterator> { protected: using iterator = ::akantu::types::details::vector_iterator; public: using iterator_category = typename iterator::iterator_category; using value_type = typename iterator::value_type; using difference_type = typename iterator::difference_type; using pointer = typename iterator::pointer; using reference = typename iterator::reference; }; + +template +struct iterator_traits<::akantu::types::details::column_iterator> { +protected: + using iterator = ::akantu::types::details::column_iterator; + +public: + using iterator_category = typename iterator::iterator_category; + using value_type = typename iterator::value_type; + using difference_type = typename iterator::difference_type; + using pointer = typename iterator::pointer; + using reference = typename iterator::reference; +}; + } // namespace std #endif /* __AKANTU_AKA_TYPES_HH__ */ diff --git a/test/test_common/test_tensors.cc b/test/test_common/test_tensors.cc index 12427290c..01acc27e6 100644 --- a/test/test_common/test_tensors.cc +++ b/test/test_common/test_tensors.cc @@ -1,578 +1,594 @@ /** * @file test_tensors.cc * * @author Nicolas Richart * * @date creation: Tue Nov 14 2017 * @date last modification: Mon Jan 22 2018 * * @brief test the tensors types * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_types.hh" +#include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { /* -------------------------------------------------------------------------- */ class TensorConstructorFixture : public ::testing::Test { public: void SetUp() override { for (auto & r : reference) { r = rand(); // google-test seeds srand() } } void TearDown() override {} template void compareToRef(const V & v) { for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i], v.storage()[i]); } } protected: const int size_{24}; const std::array mat_size{{4, 6}}; // const std::array tens3_size{{4, 2, 3}}; std::array reference; }; /* -------------------------------------------------------------------------- */ class TensorFixture : public TensorConstructorFixture { public: TensorFixture() : vref(reference.data(), size_), mref(reference.data(), mat_size[0], mat_size[1]) {} protected: Vector vref; Matrix mref; }; /* -------------------------------------------------------------------------- */ // Vector ---------------------------------------------------------------------- TEST_F(TensorConstructorFixture, VectorDefaultConstruct) { Vector v; EXPECT_EQ(0, v.size()); EXPECT_EQ(nullptr, v.storage()); EXPECT_EQ(false, v.isWrapped()); } TEST_F(TensorConstructorFixture, VectorConstruct1) { double r = rand(); Vector v(size_, r); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(r, v(i)); EXPECT_DOUBLE_EQ(r, v[i]); } } TEST_F(TensorConstructorFixture, VectorConstructWrapped) { Vector v(reference.data(), size_); EXPECT_EQ(size_, v.size()); EXPECT_EQ(true, v.isWrapped()); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i], v(i)); EXPECT_DOUBLE_EQ(reference[i], v[i]); } } TEST_F(TensorConstructorFixture, VectorConstructInitializer) { Vector v{0., 1., 2., 3., 4., 5.}; EXPECT_EQ(6, v.size()); EXPECT_EQ(false, v.isWrapped()); for (int i = 0; i < 6; ++i) { EXPECT_DOUBLE_EQ(i, v(i)); } } TEST_F(TensorConstructorFixture, VectorConstructCopy1) { Vector vref(reference.data(), reference.size()); Vector v(vref); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); compareToRef(v); } TEST_F(TensorConstructorFixture, VectorConstructCopy2) { Vector vref(reference.data(), reference.size()); Vector v(vref, false); EXPECT_EQ(size_, v.size()); EXPECT_EQ(true, v.isWrapped()); compareToRef(v); } TEST_F(TensorConstructorFixture, VectorConstructProxy1) { VectorProxy vref(reference.data(), reference.size()); EXPECT_EQ(size_, vref.size()); compareToRef(vref); Vector v(vref); EXPECT_EQ(size_, v.size()); EXPECT_EQ(true, v.isWrapped()); compareToRef(v); } TEST_F(TensorConstructorFixture, VectorConstructProxy2) { Vector vref(reference.data(), reference.size()); VectorProxy v(vref); EXPECT_EQ(size_, v.size()); compareToRef(v); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorEqual) { Vector v; v = vref; compareToRef(v); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); } TEST_F(TensorFixture, VectorEqualProxy) { VectorProxy vref_proxy(vref); Vector v; v = vref; compareToRef(v); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); } TEST_F(TensorFixture, VectorEqualProxy2) { Vector v_store(size_, 0.); VectorProxy v(v_store); v = vref; compareToRef(v); compareToRef(v_store); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorSet) { Vector v(vref); compareToRef(v); double r = rand(); v.set(r); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(r, v[i]); } TEST_F(TensorFixture, VectorClear) { Vector v(vref); compareToRef(v); v.clear(); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0, v[i]); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorDivide) { Vector v; double r = rand(); v = vref / r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, v[i]); } TEST_F(TensorFixture, VectorMultiply1) { Vector v; double r = rand(); v = vref * r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } TEST_F(TensorFixture, VectorMultiply2) { Vector v; double r = rand(); v = r * vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } TEST_F(TensorFixture, VectorAddition) { Vector v; v = vref + vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., v[i]); } TEST_F(TensorFixture, VectorSubstract) { Vector v; v = vref - vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., v[i]); } TEST_F(TensorFixture, VectorDivideEqual) { Vector v(vref); double r = rand(); v /= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, v[i]); } TEST_F(TensorFixture, VectorMultiplyEqual1) { Vector v(vref); double r = rand(); v *= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } TEST_F(TensorFixture, VectorMultiplyEqual2) { Vector v(vref); v *= v; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * reference[i], v[i]); } TEST_F(TensorFixture, VectorAdditionEqual) { Vector v(vref); v += vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., v[i]); } TEST_F(TensorFixture, VectorSubstractEqual) { Vector v(vref); v -= vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., v[i]); } /* -------------------------------------------------------------------------- */ // Matrix ---------------------------------------------------------------------- TEST_F(TensorConstructorFixture, MatrixDefaultConstruct) { Matrix m; EXPECT_EQ(0, m.size()); EXPECT_EQ(0, m.rows()); EXPECT_EQ(0, m.cols()); EXPECT_EQ(nullptr, m.storage()); EXPECT_EQ(false, m.isWrapped()); } TEST_F(TensorConstructorFixture, MatrixConstruct1) { double r = rand(); Matrix m(mat_size[0], mat_size[1], r); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); for (int i = 0; i < mat_size[0]; ++i) { for (int j = 0; j < mat_size[1]; ++j) { EXPECT_EQ(r, m(i, j)); EXPECT_EQ(r, m[i + j * mat_size[0]]); } } } TEST_F(TensorConstructorFixture, MatrixConstructWrapped) { Matrix m(reference.data(), mat_size[0], mat_size[1]); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(true, m.isWrapped()); for (int i = 0; i < mat_size[0]; ++i) { for (int j = 0; j < mat_size[1]; ++j) { EXPECT_DOUBLE_EQ(reference[i + j * mat_size[0]], m(i, j)); } } compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructInitializer) { Matrix m{{0., 1., 2.}, {3., 4., 5.}}; EXPECT_EQ(6, m.size()); EXPECT_EQ(2, m.rows()); EXPECT_EQ(3, m.cols()); EXPECT_EQ(false, m.isWrapped()); int c = 0; for (int i = 0; i < 2; ++i) { for (int j = 0; j < 3; ++j, ++c) { EXPECT_DOUBLE_EQ(c, m(i, j)); } } } TEST_F(TensorConstructorFixture, MatrixConstructCopy1) { Matrix mref(reference.data(), mat_size[0], mat_size[1]); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructCopy2) { Matrix mref(reference.data(), mat_size[0], mat_size[1]); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructProxy1) { MatrixProxy mref(reference.data(), mat_size[0], mat_size[1]); EXPECT_EQ(size_, mref.size()); EXPECT_EQ(mat_size[0], mref.size(0)); EXPECT_EQ(mat_size[1], mref.size(1)); compareToRef(mref); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(true, m.isWrapped()); compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructProxy2) { Matrix mref(reference.data(), mat_size[0], mat_size[1]); MatrixProxy m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.size(0)); EXPECT_EQ(mat_size[1], m.size(1)); compareToRef(m); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixEqual) { Matrix m; m = mref; compareToRef(m); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); } TEST_F(TensorFixture, MatrixEqualProxy1) { MatrixProxy mref_proxy(mref); Matrix m; m = mref; compareToRef(m); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); } TEST_F(TensorFixture, MatrixEqualProxy2) { Matrix m_store(mat_size[0], mat_size[1], 0.); MatrixProxy m(m_store); m = mref; compareToRef(m); compareToRef(m_store); } TEST_F(TensorFixture, MatrixEqualSlice) { Matrix m(mat_size[0], mat_size[1], 0.); for (unsigned int i = 0; i < m.cols(); ++i) m(i) = Vector(mref(i)); compareToRef(m); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixSet) { Matrix m(mref); compareToRef(m); double r = rand(); m.set(r); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(r, m[i]); } TEST_F(TensorFixture, MatrixClear) { Matrix m(mref); compareToRef(m); m.clear(); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0, m[i]); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixDivide) { Matrix m; double r = rand(); m = mref / r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, m[i]); } TEST_F(TensorFixture, MatrixMultiply1) { Matrix m; double r = rand(); m = mref * r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, m[i]); } TEST_F(TensorFixture, MatrixMultiply2) { Matrix m; double r = rand(); m = r * mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, m[i]); } TEST_F(TensorFixture, MatrixAddition) { Matrix m; m = mref + mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., m[i]); } TEST_F(TensorFixture, MatrixSubstract) { Matrix m; m = mref - mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., m[i]); } TEST_F(TensorFixture, MatrixDivideEqual) { Matrix m(mref); double r = rand(); m /= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, m[i]); } TEST_F(TensorFixture, MatrixMultiplyEqual1) { Matrix m(mref); double r = rand(); m *= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, m[i]); } TEST_F(TensorFixture, MatrixAdditionEqual) { Matrix m(mref); m += mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., m[i]); } TEST_F(TensorFixture, MatrixSubstractEqual) { Matrix m(mref); m -= mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., m[i]); } TEST_F(TensorFixture, MatrixIterator) { Matrix m(mref); UInt col_count = 0; for (auto && col : m) { Vector col_hand(m.storage() + col_count * m.rows(), m.rows()); Vector col_wrap(col); auto comp = (col_wrap - col_hand).norm(); EXPECT_DOUBLE_EQ(0., comp); ++col_count; } } +TEST_F(TensorFixture, MatrixIteratorZip) { + Matrix m1(mref); + Matrix m2(mref); + + UInt col_count = 0; + for (auto && col : zip(m1, m2)) { + Vector col1(std::get<0>(col)); + Vector col2(std::get<1>(col)); + + auto comp = (col1 - col2).norm(); + EXPECT_DOUBLE_EQ(0., comp); + ++col_count; + } +} + #if defined(AKANTU_USE_LAPACK) TEST_F(TensorFixture, MatrixEigs) { Matrix m{{0, 1, 0, 0}, {1., 0, 0, 0}, {0, 1, 0, 1}, {0, 0, 4, 0}}; Matrix eig_vects(4, 4); Vector eigs(4); m.eig(eigs, eig_vects); Vector eigs_ref{2, 1., -1., -2}; auto lambda_v = m * eig_vects; for (int i = 0; i < 4; ++i) { EXPECT_NEAR(eigs_ref(i), eigs(i), 1e-14); for (int j = 0; j < 4; ++j) { EXPECT_NEAR(lambda_v(i)(j), eigs(i) * eig_vects(i)(j), 1e-14); } } } #endif /* -------------------------------------------------------------------------- */ } // namespace