diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh
index d2b6f4e38..16a0ac321 100644
--- a/src/common/aka_array.hh
+++ b/src/common/aka_array.hh
@@ -1,363 +1,363 @@
/**
* @file aka_array.hh
*
* @author Till Junge <till.junge@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Jan 22 2016
*
* @brief Array container for Akantu
* This container differs from the std::vector from the fact it as 2 dimensions
* a main dimension and the size stored per entries
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_VECTOR_HH__
#define __AKANTU_VECTOR_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
#include <typeinfo>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace akantu {
/// class that afford to store vectors in static memory
class ArrayBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
explicit ArrayBase(ID id = "");
virtual ~ArrayBase();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// get the amount of space allocated in bytes
inline UInt getMemorySize() const;
/// set the size to zero without freeing the allocated space
inline void empty();
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// Get the real size allocated in memory
AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt);
/// Get the Size of the Array
- UInt getsize() const __attribute__((deprecated)) { return size_; }
+ UInt getSize() const __attribute__((deprecated)) { return size_; }
UInt size() const { return size_; }
/// Get the number of components
AKANTU_GET_MACRO(NbComponent, nb_component, UInt);
/// Get the name of th array
AKANTU_GET_MACRO(ID, id, const ID &);
/// Set the name of th array
AKANTU_SET_MACRO(ID, id, const ID &);
- // AKANTU_GET_MACRO(Tag, tag, const std::string &);
- // AKANTU_SET_MACRO(Tag, tag, const std::string &);
-
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// id of the vector
ID id;
/// the size allocated
UInt allocated_size{0};
/// the size used
UInt size_{0};
/// number of components
UInt nb_component{1};
/// size of the stored type
UInt size_of_type{0};
};
/* -------------------------------------------------------------------------- */
namespace {
template <std::size_t dim, typename T> struct IteratorHelper {};
template <typename T> struct IteratorHelper<0, T> { using type = T; };
template <typename T> struct IteratorHelper<1, T> { using type = Vector<T>; };
template <typename T> struct IteratorHelper<2, T> { using type = Matrix<T>; };
template <typename T> struct IteratorHelper<3, T> {
using type = Tensor3<T>;
};
template <std::size_t dim, typename T>
using IteratorHelper_t = typename IteratorHelper<dim, T>::type;
} // namespace
/* -------------------------------------------------------------------------- */
template <typename T, bool is_scal> class Array : public ArrayBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using value_type = T;
using reference = value_type &;
using pointer_type = value_type *;
using const_reference = const value_type &;
/// Allocation of a new vector
explicit inline Array(UInt size = 0, UInt nb_component = 1,
const ID & id = "");
/// Allocation of a new vector with a default value
Array(UInt size, UInt nb_component, const value_type def_values[],
const ID & id = "");
/// Allocation of a new vector with a default value
Array(UInt size, UInt nb_component, const_reference value,
const ID & id = "");
/// Copy constructor (deep copy if deep=true)
Array(const Array<value_type, is_scal> & vect, bool deep = true,
const ID & id = "");
#ifndef SWIG
/// Copy constructor (deep copy)
explicit Array(const std::vector<value_type> & vect);
#endif
inline ~Array() override;
Array & operator=(const Array & a) {
/// this is to let STL allocate and copy arrays in the case of
/// std::vector::resize
AKANTU_DEBUG_ASSERT(this->size == 0, "Cannot copy akantu::Array");
return const_cast<Array &>(a);
}
/* ------------------------------------------------------------------------ */
/* Iterator */
/* ------------------------------------------------------------------------ */
/// \todo protected: does not compile with intel check why
public:
- template <class R, class IR = R, bool issame = std::is_same<IR, T>::value>
+ template <class R, class it, class IR = R, bool is_tensor_ = is_tensor<R>{}>
class iterator_internal;
public:
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
template <typename R = T> class const_iterator;
template <typename R = T> class iterator;
/* ------------------------------------------------------------------------ */
/// iterator for Array of nb_component = 1
using scalar_iterator = iterator<T>;
/// const_iterator for Array of nb_component = 1
using const_scalar_iterator = const_iterator<T>;
/// iterator returning Vectors of size n on entries of Array with
/// nb_component = n
using vector_iterator = iterator<Vector<T>>;
/// const_iterator returning Vectors of n size on entries of Array with
/// nb_component = n
using const_vector_iterator = const_iterator<Vector<T>>;
/// iterator returning Matrices of size (m, n) on entries of Array with
/// nb_component = m*n
using matrix_iterator = iterator<Matrix<T>>;
/// const iterator returning Matrices of size (m, n) on entries of Array with
/// nb_component = m*n
using const_matrix_iterator = const_iterator<Matrix<T>>;
/// iterator returning Tensor3 of size (m, n, k) on entries of Array with
/// nb_component = m*n*k
using tensor3_iterator = iterator<Tensor3<T>>;
/// const iterator returning Tensor3 of size (m, n, k) on entries of Array
/// with nb_component = m*n*k
using const_tensor3_iterator = const_iterator<Tensor3<T>>;
/* ------------------------------------------------------------------------ */
- template <typename... Ns> inline decltype(auto) begin(Ns... n);
- template <typename... Ns> inline decltype(auto) end(Ns... n);
+ template <typename... Ns> inline decltype(auto) begin(Ns&&... n);
+ template <typename... Ns> inline decltype(auto) end(Ns&&... n);
- template <typename... Ns> inline decltype(auto) begin(Ns... n) const;
- template <typename... Ns> inline decltype(auto) end(Ns... n) const;
+ template <typename... Ns> inline decltype(auto) begin(Ns&&... n) const;
+ template <typename... Ns> inline decltype(auto) end(Ns&&... n) const;
- template <typename... Ns>
- inline decltype(auto) begin_reinterpret(Ns... n);
- template <typename... Ns>
- inline decltype(auto) end_reinterpret(Ns... n);
+ template <typename... Ns> inline decltype(auto) begin_reinterpret(Ns&&... n);
+ template <typename... Ns> inline decltype(auto) end_reinterpret(Ns&&... n);
template <typename... Ns>
- inline decltype(auto) begin_reinterpret(Ns... n) const;
+ inline decltype(auto) begin_reinterpret(Ns&&... n) const;
template <typename... Ns>
- inline decltype(auto) end_reinterpret(Ns... n) const;
+ inline decltype(auto) end_reinterpret(Ns&&... n) const;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// append a tuple of size nb_component containing value
inline void push_back(const_reference value);
/// append a vector
// inline void push_back(const value_type new_elem[]);
/// append a Vector or a Matrix
- template <template <typename> class C>
+ template <template <typename> class C,
+ typename = std::enable_if_t<is_tensor<C<T>>{}>>
inline void push_back(const C<T> & new_elem);
/// append the value of the iterator
template <typename Ret> inline void push_back(const iterator<Ret> & it);
/// erase the value at position i
inline void erase(UInt i);
/// ask Nico, clarify
template <typename R> inline iterator<R> erase(const iterator<R> & it);
/// changes the allocated size but not the size
virtual void reserve(UInt size);
/// change the size of the Array
virtual void resize(UInt size);
/// change the size of the Array and initialize the values
virtual void resize(UInt size, const T & val);
/// change the number of components by interlacing data
/// @param multiplicator number of interlaced components add
/// @param block_size blocks of data in the array
/// Examaple for block_size = 2, multiplicator = 2
/// array = oo oo oo -> new array = oo nn nn oo nn nn oo nn nn
void extendComponentsInterlaced(UInt multiplicator, UInt stride);
/// search elem in the vector, return the position of the first occurrence or
/// -1 if not found
UInt find(const_reference elem) const;
/// @see Array::find(const_reference elem) const
UInt find(T elem[]) const;
/// @see Array::find(const_reference elem) const
- template <template <typename> class C> inline UInt find(const C<T> & elem);
+ template <template <typename> class C,
+ typename = std::enable_if_t<is_tensor<C<T>>{}>>
+ inline UInt find(const C<T> & elem);
/// set all entries of the array to 0
inline void clear() { std::fill_n(values, size_ * nb_component, T()); }
/// set all entries of the array to the value t
/// @param t value to fill the array with
inline void set(T t) { std::fill_n(values, size_ * nb_component, t); }
/// set all tuples of the array to a given vector or matrix
/// @param vm Matrix or Vector to fill the array with
- template <template <typename> class C> inline void set(const C<T> & vm);
+ template <template <typename> class C,
+ typename = std::enable_if_t<is_tensor<C<T>>{}>>
+ inline void set(const C<T> & vm);
/// Append the content of the other array to the current one
void append(const Array<T> & other);
/// copy another Array in the current Array, the no_sanity_check allows you to
/// force the copy in cases where you know what you do with two non matching
/// Arrays in terms of n
void copy(const Array<T, is_scal> & other, bool no_sanity_check = false);
/// give the address of the memory allocated for this vector
T * storage() const { return values; };
/// function to print the containt of the class
void printself(std::ostream & stream, int indent = 0) const override;
protected:
/// perform the allocation for the constructors
void allocate(UInt size, UInt nb_component = 1);
/// resize initializing with uninitialized_fill if fill is set
void resizeUnitialized(UInt new_size, bool fill, const T & val = T());
/* ------------------------------------------------------------------------ */
/* Operators */
/* ------------------------------------------------------------------------ */
public:
/// substraction entry-wise
Array<T, is_scal> & operator-=(const Array<T, is_scal> & other);
/// addition entry-wise
Array<T, is_scal> & operator+=(const Array<T, is_scal> & other);
/// multiply evry entry by alpha
Array<T, is_scal> & operator*=(const T & alpha);
/// check if the array are identical entry-wise
bool operator==(const Array<T, is_scal> & other) const;
/// @see Array::operator==(const Array<T, is_scal> & other) const
bool operator!=(const Array<T, is_scal> & other) const;
/// return a reference to the j-th entry of the i-th tuple
inline reference operator()(UInt i, UInt j = 0);
/// return a const reference to the j-th entry of the i-th tuple
inline const_reference operator()(UInt i, UInt j = 0) const;
/// return a reference to the ith component of the 1D array
inline reference operator[](UInt i);
/// return a const reference to the ith component of the 1D array
inline const_reference operator[](UInt i) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// array of values
T * values; // /!\ very dangerous
};
/* -------------------------------------------------------------------------- */
/* Inline Functions Array<T, is_scal> */
/* -------------------------------------------------------------------------- */
template <typename T, bool is_scal>
inline std::ostream & operator<<(std::ostream & stream,
const Array<T, is_scal> & _this) {
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
/* Inline Functions ArrayBase */
/* -------------------------------------------------------------------------- */
inline std::ostream & operator<<(std::ostream & stream,
const ArrayBase & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "aka_array_tmpl.hh"
#include "aka_types.hh"
#endif /* __AKANTU_VECTOR_HH__ */
diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh
index 41c10a307..eb9b21b7e 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1251 +1,1252 @@
/**
* @file aka_array_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Jul 15 2010
* @date last modification: Fri Jan 22 2016
*
* @brief Inline functions of the classes Array<T> and ArrayBase
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* Inline Functions Array<T> */
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
/* -------------------------------------------------------------------------- */
#include <memory>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__
#define __AKANTU_AKA_ARRAY_TMPL_HH__
namespace akantu {
+namespace debug {
+ struct ArrayException : public Exception {};
+} // namespace debug
+
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline T & Array<T, is_scal>::operator()(UInt i, UInt j) {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
"The value at position ["
<< i << "," << j << "] is out of range in array \""
<< id << "\"");
return values[i * nb_component + j];
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline const T & Array<T, is_scal>::operator()(UInt i, UInt j) const {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
"The value at position ["
<< i << "," << j << "] is out of range in array \""
<< id << "\"");
return values[i * nb_component + j];
}
template <class T, bool is_scal>
inline T & Array<T, is_scal>::operator[](UInt i) {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
"The value at position ["
<< i << "] is out of range in array \"" << id
<< "\"");
return values[i];
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline const T & Array<T, is_scal>::operator[](UInt i) const {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
"The value at position ["
<< i << "] is out of range in array \"" << id
<< "\"");
return values[i];
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array with the value value for all components
* @param value the new last tuple or the array will contain nb_component copies
* of value
*/
template <class T, bool is_scal>
inline void Array<T, is_scal>::push_back(const T & value) {
resizeUnitialized(size_ + 1, true, value);
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array
* @param new_elem a C-array containing the values to be copied to the end of
* the array */
// template <class T, bool is_scal>
// inline void Array<T, is_scal>::push_back(const T new_elem[]) {
// UInt pos = size_;
// resizeUnitialized(size_ + 1, false);
// T * tmp = values + nb_component * pos;
// std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
// }
/* -------------------------------------------------------------------------- */
/**
* append a matrix or a vector to the array
* @param new_elem a reference to a Matrix<T> or Vector<T> */
template <class T, bool is_scal>
-template <template <typename> class C>
+template <template <typename> class C, typename>
inline void Array<T, is_scal>::push_back(const C<T> & new_elem) {
AKANTU_DEBUG_ASSERT(
nb_component == new_elem.size(),
"The vector("
<< new_elem.size()
<< ") as not a size compatible with the Array (nb_component="
<< nb_component << ").");
UInt pos = size_;
resizeUnitialized(size_ + 1, false);
T * tmp = values + nb_component * pos;
std::uninitialized_copy(new_elem.storage(), new_elem.storage() + nb_component,
tmp);
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array
* @param it an iterator to the tuple to be copied to the end of the array */
template <class T, bool is_scal>
template <class Ret>
inline void
Array<T, is_scal>::push_back(const Array<T, is_scal>::iterator<Ret> & it) {
UInt pos = size_;
resizeUnitialized(size_ + 1, false);
T * tmp = values + nb_component * pos;
T * new_elem = it.data();
std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
}
/* -------------------------------------------------------------------------- */
/**
* erase an element. If the erased element is not the last of the array, the
* last element is moved into the hole in order to maintain contiguity. This
* may invalidate existing iterators (For instance an iterator obtained by
* Array::end() is no longer correct) and will change the order of the
* elements.
* @param i index of element to erase
*/
template <class T, bool is_scal> inline void Array<T, is_scal>::erase(UInt i) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT((size_ > 0), "The array is empty");
AKANTU_DEBUG_ASSERT((i < size_), "The element at position ["
<< i << "] is out of range (" << i
<< ">=" << size_ << ")");
if (i != (size_ - 1)) {
for (UInt j = 0; j < nb_component; ++j) {
values[i * nb_component + j] = values[(size_ - 1) * nb_component + j];
}
}
resize(size_ - 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Subtract another array entry by entry from this array in place. Both arrays
* must
* have the same size and nb_component. If the arrays have different shapes,
* code compiled in debug mode will throw an expeption and optimised code
* will behave in an unpredicted manner
* @param other array to subtract from this
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::
operator-=(const Array<T, is_scal> & vect) {
AKANTU_DEBUG_ASSERT((size_ == vect.size_) &&
(nb_component == vect.nb_component),
"The too array don't have the same sizes");
T * a = values;
T * b = vect.storage();
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a -= *b;
++a;
++b;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Add another array entry by entry to this array in place. Both arrays must
* have the same size and nb_component. If the arrays have different shapes,
* code compiled in debug mode will throw an expeption and optimised code
* will behave in an unpredicted manner
* @param other array to add to this
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::
operator+=(const Array<T, is_scal> & vect) {
AKANTU_DEBUG_ASSERT((size_ == vect.size) &&
(nb_component == vect.nb_component),
"The too array don't have the same sizes");
T * a = values;
T * b = vect.storage();
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a++ += *b++;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Multiply all entries of this array by a scalar in place
* @param alpha scalar multiplicant
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::operator*=(const T & alpha) {
T * a = values;
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a++ *= alpha;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Compare this array element by element to another.
* @param other array to compare to
* @return true it all element are equal and arrays have the same shape, else
* false
*/
template <class T, bool is_scal>
bool Array<T, is_scal>::operator==(const Array<T, is_scal> & array) const {
bool equal = nb_component == array.nb_component && size_ == array.size_ &&
id == array.id;
if (!equal)
return false;
if (values == array.storage())
return true;
else
return std::equal(values, values + size_ * nb_component, array.storage());
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
bool Array<T, is_scal>::operator!=(const Array<T, is_scal> & array) const {
return !operator==(array);
}
/* -------------------------------------------------------------------------- */
/**
* set all tuples of the array to a given vector or matrix
* @param vm Matrix or Vector to fill the array with
*/
template <class T, bool is_scal>
-template <template <typename> class C>
+template <template <typename> class C, typename>
inline void Array<T, is_scal>::set(const C<T> & vm) {
AKANTU_DEBUG_ASSERT(
nb_component == vm.size(),
"The size of the object does not match the number of components");
for (T * it = values; it < values + nb_component * size_;
it += nb_component) {
std::copy(vm.storage(), vm.storage() + nb_component, it);
}
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::append(const Array<T> & other) {
AKANTU_DEBUG_ASSERT(
nb_component == other.nb_component,
"Cannot append an array with a different number of component");
UInt old_size = this->size_;
this->resizeUnitialized(this->size_ + other.size(), false);
T * tmp = values + nb_component * old_size;
std::uninitialized_copy(other.storage(),
other.storage() + other.size() * nb_component, tmp);
}
/* -------------------------------------------------------------------------- */
/* Functions Array<T, is_scal> */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const ID & id)
: ArrayBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
if (!is_scal) {
T val = T();
std::uninitialized_fill(values, values + size * nb_component, val);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const T def_values[],
const ID & id)
: ArrayBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
T * tmp = values;
for (UInt i = 0; i < size; ++i) {
tmp = values + nb_component * i;
std::uninitialized_copy(def_values, def_values + nb_component, tmp);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const T & value,
const ID & id)
: ArrayBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
std::uninitialized_fill_n(values, size * nb_component, value);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(const Array<T, is_scal> & vect, bool deep,
const ID & id)
: ArrayBase(vect) {
AKANTU_DEBUG_IN();
this->id = (id == "") ? vect.id : id;
if (deep) {
allocate(vect.size_, vect.nb_component);
T * tmp = values;
std::uninitialized_copy(vect.storage(),
vect.storage() + size_ * nb_component, tmp);
} else {
this->values = vect.storage();
this->size_ = vect.size_;
this->nb_component = vect.nb_component;
this->allocated_size = vect.allocated_size;
this->size_of_type = vect.size_of_type;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
#ifndef SWIG
template <class T, bool is_scal>
Array<T, is_scal>::Array(const std::vector<T> & vect) {
AKANTU_DEBUG_IN();
this->id = "";
allocate(vect.size(), 1);
T * tmp = values;
std::uninitialized_copy(&(vect[0]), &(vect[size_ - 1]), tmp);
AKANTU_DEBUG_OUT();
}
#endif
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> Array<T, is_scal>::~Array() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG(dblAccessory,
"Freeing " << printMemorySize<T>(allocated_size * nb_component)
<< " (" << id << ")");
if (values) {
if (!is_scal)
for (UInt i = 0; i < size_ * nb_component; ++i) {
T * obj = values + i;
obj->~T();
}
free(values);
}
size_ = allocated_size = 0;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* perform the allocation for the constructors
* @param size is the size of the array
* @param nb_component is the number of component of the array
*/
template <class T, bool is_scal>
void Array<T, is_scal>::allocate(UInt size, UInt nb_component) {
AKANTU_DEBUG_IN();
if (size == 0) {
values = nullptr;
} else {
values = static_cast<T *>(malloc(nb_component * size * sizeof(T)));
AKANTU_DEBUG_ASSERT(values != nullptr,
"Cannot allocate "
<< printMemorySize<T>(size * nb_component) << " ("
<< id << ")");
}
if (values == NULL) {
this->size_ = this->allocated_size = 0;
} else {
AKANTU_DEBUG(dblAccessory, "Allocated "
<< printMemorySize<T>(size * nb_component)
<< " (" << id << ")");
this->size_ = this->allocated_size = size;
}
this->size_of_type = sizeof(T);
this->nb_component = nb_component;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::reserve(UInt new_size) {
UInt tmp_size = this->size_;
resizeUnitialized(new_size, false);
this->size_ = tmp_size;
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed. If the
* size increases, the new tuples are filled with zeros
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal> void Array<T, is_scal>::resize(UInt new_size) {
resizeUnitialized(new_size, !is_scal);
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed. If the
* size increases, the new tuples are filled with zeros
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal>
void Array<T, is_scal>::resize(UInt new_size, const T & val) {
this->resizeUnitialized(new_size, true, val);
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed.
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal>
void Array<T, is_scal>::resizeUnitialized(UInt new_size, bool fill,
const T & val) {
// AKANTU_DEBUG_IN();
// free some memory
if (new_size <= allocated_size) {
if (!is_scal) {
T * old_values = values;
if (new_size < size_) {
for (UInt i = new_size * nb_component; i < size_ * nb_component; ++i) {
T * obj = old_values + i;
obj->~T();
}
}
}
if (allocated_size - new_size > AKANTU_MIN_ALLOCATION) {
AKANTU_DEBUG(dblAccessory,
"Freeing " << printMemorySize<T>((allocated_size - size_) *
nb_component)
<< " (" << id << ")");
// Normally there are no allocation problem when reducing an array
if (new_size == 0) {
free(values);
- values = NULL;
+ values = nullptr;
} else {
auto * tmp_ptr = static_cast<T *>(
realloc(values, new_size * nb_component * sizeof(T)));
- if (tmp_ptr == NULL) {
+ if (tmp_ptr == nullptr) {
AKANTU_EXCEPTION("Cannot free data ("
<< id << ")"
<< " [current allocated size : " << allocated_size
<< " | "
<< "requested size : " << new_size << "]");
}
values = tmp_ptr;
}
allocated_size = new_size;
}
} else {
// allocate more memory
UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION)
? allocated_size + AKANTU_MIN_ALLOCATION
: new_size;
auto * tmp_ptr = static_cast<T *>(
realloc(values, size_to_alloc * nb_component * sizeof(T)));
AKANTU_DEBUG_ASSERT(
tmp_ptr != NULL,
"Cannot allocate " << printMemorySize<T>(size_to_alloc * nb_component));
if (tmp_ptr == NULL) {
AKANTU_DEBUG_ERROR("Cannot allocate more data ("
<< id << ")"
<< " [current allocated size : " << allocated_size
<< " | "
<< "requested size : " << new_size << "]");
}
AKANTU_DEBUG(dblAccessory,
"Allocating " << printMemorySize<T>(
(size_to_alloc - allocated_size) * nb_component));
allocated_size = size_to_alloc;
values = tmp_ptr;
}
if (fill && this->size_ < new_size) {
std::uninitialized_fill(values + size_ * nb_component,
values + new_size * nb_component, val);
}
size_ = new_size;
// AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* change the number of components by interlacing data
* @param multiplicator number of interlaced components add
* @param block_size blocks of data in the array
* Examaple for block_size = 2, multiplicator = 2
* array = oo oo oo -> new array = oo nn nn oo nn nn oo nn nn */
template <class T, bool is_scal>
void Array<T, is_scal>::extendComponentsInterlaced(UInt multiplicator,
UInt block_size) {
AKANTU_DEBUG_IN();
if (multiplicator == 1)
return;
AKANTU_DEBUG_ASSERT(multiplicator > 1, "invalid multiplicator");
AKANTU_DEBUG_ASSERT(nb_component % block_size == 0,
"stride must divide actual number of components");
values = static_cast<T *>(
realloc(values, nb_component * multiplicator * size_ * sizeof(T)));
UInt new_component = nb_component / block_size * multiplicator;
for (UInt i = 0, k = size_ - 1; i < size_; ++i, --k) {
for (UInt j = 0; j < new_component; ++j) {
UInt m = new_component - j - 1;
UInt n = m / multiplicator;
for (UInt l = 0, p = block_size - 1; l < block_size; ++l, --p) {
values[k * nb_component * multiplicator + m * block_size + p] =
values[k * nb_component + n * block_size + p];
}
}
}
nb_component = nb_component * multiplicator;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* search elem in the array, return the position of the first occurrence or
* -1 if not found
* @param elem the element to look for
* @return index of the first occurrence of elem or -1 if elem is not present
*/
template <class T, bool is_scal>
UInt Array<T, is_scal>::find(const T & elem) const {
AKANTU_DEBUG_IN();
auto begin = this->begin();
auto end = this->end();
auto it = std::find(begin, end, elem);
AKANTU_DEBUG_OUT();
return (it != end) ? it - begin : UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> UInt Array<T, is_scal>::find(T elem[]) const {
AKANTU_DEBUG_IN();
T * it = values;
UInt i = 0;
for (; i < size_; ++i) {
if (*it == elem[0]) {
T * cit = it;
UInt c = 0;
for (; (c < nb_component) && (*cit == elem[c]); ++c, ++cit)
;
if (c == nb_component) {
AKANTU_DEBUG_OUT();
return i;
}
}
it += nb_component;
}
return UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
-template <template <typename> class C>
+template <template <typename> class C, typename>
inline UInt Array<T, is_scal>::find(const C<T> & elem) {
AKANTU_DEBUG_ASSERT(elem.size() == nb_component,
"Cannot find an element with a wrong size ("
<< elem.size() << ") != " << nb_component);
return this->find(elem.storage());
}
/* -------------------------------------------------------------------------- */
/**
* copy the content of another array. This overwrites the current content.
* @param other Array to copy into this array. It has to have the same
* nb_component as this. If compiled in debug mode, an incorrect other will
* result in an exception being thrown. Optimised code may result in
* unpredicted behaviour.
*/
template <class T, bool is_scal>
void Array<T, is_scal>::copy(const Array<T, is_scal> & vect,
bool no_sanity_check) {
AKANTU_DEBUG_IN();
if (!no_sanity_check)
if (vect.nb_component != nb_component)
AKANTU_DEBUG_ERROR(
"The two arrays do not have the same number of components");
resize((vect.size_ * vect.nb_component) / nb_component);
T * tmp = values;
std::uninitialized_copy(vect.storage(), vect.storage() + size_ * nb_component,
tmp);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <bool is_scal> class ArrayPrintHelper {
public:
template <typename T>
static void print_content(const Array<T> & vect, std::ostream & stream,
int indent) {
if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << " + values : {";
for (UInt i = 0; i < vect.size(); ++i) {
stream << "{";
for (UInt j = 0; j < vect.getNbComponent(); ++j) {
stream << vect(i, j);
if (j != vect.getNbComponent() - 1)
stream << ", ";
}
stream << "}";
if (i != vect.size() - 1)
stream << ", ";
}
stream << "}" << std::endl;
}
}
};
template <> class ArrayPrintHelper<false> {
public:
template <typename T>
static void print_content(__attribute__((unused)) const Array<T> & vect,
__attribute__((unused)) std::ostream & stream,
__attribute__((unused)) int indent) {}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
std::streamsize prec = stream.precision();
std::ios_base::fmtflags ff = stream.flags();
stream.setf(std::ios_base::showbase);
stream.precision(2);
stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> ["
<< std::endl;
stream << space << " + id : " << this->id << std::endl;
stream << space << " + size : " << this->size_ << std::endl;
stream << space << " + nb_component : " << this->nb_component << std::endl;
stream << space << " + allocated size : " << this->allocated_size
<< std::endl;
stream << space << " + memory size : "
<< printMemorySize<T>(allocated_size * nb_component) << std::endl;
if (!AKANTU_DEBUG_LEVEL_IS_TEST())
stream << space << " + address : " << std::hex << this->values
<< std::dec << std::endl;
stream.precision(prec);
stream.flags(ff);
ArrayPrintHelper<is_scal>::print_content(*this, stream, indent);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
/* Inline Functions ArrayBase */
/* -------------------------------------------------------------------------- */
inline UInt ArrayBase::getMemorySize() const {
return allocated_size * nb_component * size_of_type;
}
inline void ArrayBase::empty() { size_ = 0; }
/* -------------------------------------------------------------------------- */
/* Iterators */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
-template <class R, class IR, bool is_r_scal>
+template <class R, class daughter, class IR, bool is_tensor>
class Array<T, is_scal>::iterator_internal {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using proxy = typename R::proxy;
using const_proxy = const typename R::proxy;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
- iterator_internal() : initial(NULL), ret(NULL), ret_ptr(NULL){};
+ iterator_internal() = default;
iterator_internal(pointer_type data, UInt _offset)
- : _offset(_offset), initial(data), ret(NULL), ret_ptr(data) {
+ : _offset(_offset), initial(data), ret(nullptr), ret_ptr(data) {
AKANTU_DEBUG_ERROR(
"The constructor should never be called it is just an ugly trick...");
}
- iterator_internal(pointer wrapped)
+ iterator_internal(std::unique_ptr<internal_value_type> && wrapped)
: _offset(wrapped->size()), initial(wrapped->storage()),
- ret(const_cast<internal_pointer>(wrapped)),
- ret_ptr(wrapped->storage()) {}
+ ret(std::move(wrapped)), ret_ptr(ret->storage()) {}
iterator_internal(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
- this->ret = new internal_value_type(*it.ret, false);
+ this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
}
- virtual ~iterator_internal() { delete ret; };
+ iterator_internal(iterator_internal && it) = default;
+
+ virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
if (this->ret)
this->ret->shallowCopy(*it.ret);
else
- this->ret = new internal_value_type(*it.ret, false);
+ this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
return *this;
}
UInt getCurrentIndex() {
return (this->ret_ptr - this->initial) / this->_offset;
};
inline reference operator*() {
ret->values = ret_ptr;
return *ret;
};
inline const_reference operator*() const {
ret->values = ret_ptr;
return *ret;
};
inline pointer operator->() {
ret->values = ret_ptr;
- return ret;
+ return ret.get();
};
- inline iterator_internal & operator++() {
+ inline daughter & operator++() {
ret_ptr += _offset;
- return *this;
+ return static_cast<daughter &>(*this);
};
- inline iterator_internal & operator--() {
+ inline daughter & operator--() {
ret_ptr -= _offset;
- return *this;
+ return static_cast<daughter &>(*this);
};
- inline iterator_internal & operator+=(const UInt n) {
+ inline daughter & operator+=(const UInt n) {
ret_ptr += _offset * n;
- return *this;
+ return static_cast<daughter &>(*this);
}
- inline iterator_internal & operator-=(const UInt n) {
+ inline daughter & operator-=(const UInt n) {
ret_ptr -= _offset * n;
- return *this;
+ return static_cast<daughter &>(*this);
}
inline proxy operator[](const UInt n) {
ret->values = ret_ptr + n * _offset;
return proxy(*ret);
}
inline const_proxy operator[](const UInt n) const {
ret->values = ret_ptr + n * _offset;
return const_proxy(*ret);
}
inline bool operator==(const iterator_internal & other) const {
return this->ret_ptr == other.ret_ptr;
}
inline bool operator!=(const iterator_internal & other) const {
return this->ret_ptr != other.ret_ptr;
}
inline bool operator<(const iterator_internal & other) const {
return this->ret_ptr < other.ret_ptr;
}
inline bool operator<=(const iterator_internal & other) const {
return this->ret_ptr <= other.ret_ptr;
}
inline bool operator>(const iterator_internal & other) const {
return this->ret_ptr > other.ret_ptr;
}
inline bool operator>=(const iterator_internal & other) const {
return this->ret_ptr >= other.ret_ptr;
}
- inline iterator_internal operator+(difference_type n) {
- iterator_internal tmp(*this);
+ inline daughter operator+(difference_type n) {
+ daughter tmp(static_cast<daughter &>(*this));
tmp += n;
return tmp;
}
- inline iterator_internal operator-(difference_type n) {
- iterator_internal tmp(*this);
+ inline daughter operator-(difference_type n) {
+ daughter tmp(static_cast<daughter &>(*this));
tmp -= n;
return tmp;
}
inline difference_type operator-(const iterator_internal & b) {
return (this->ret_ptr - b.ret_ptr) / _offset;
}
inline pointer_type data() const { return ret_ptr; }
inline difference_type offset() const { return _offset; }
protected:
UInt _offset{0};
- pointer_type initial;
- internal_pointer ret;
- pointer_type ret_ptr;
+ pointer_type initial{nullptr};
+ std::unique_ptr<internal_value_type> ret{nullptr};
+ pointer_type ret_ptr{nullptr};
};
/* -------------------------------------------------------------------------- */
/**
* Specialization for scalar types
*/
template <class T, bool is_scal>
-template <class R, class IR>
-class Array<T, is_scal>::iterator_internal<R, IR, true> {
+template <class R, class daughter, class IR>
+class Array<T, is_scal>::iterator_internal<R, daughter, IR, false> {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
- iterator_internal(pointer data = nullptr, UInt _offset = 1)
- : _offset(_offset), ret(data), initial(data){};
+ iterator_internal(pointer data = nullptr) : ret(data), initial(data){};
iterator_internal(const iterator_internal & it) = default;
iterator_internal(iterator_internal && it) = default;
virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) = default;
- UInt getCurrentIndex() {
- return (this->ret - this->initial) / this->_offset;
- };
+ UInt getCurrentIndex() { return (this->ret - this->initial); };
inline reference operator*() { return *ret; };
inline const_reference operator*() const { return *ret; };
inline pointer operator->() { return ret; };
- inline iterator_internal & operator++() {
+ inline daughter & operator++() {
++ret;
- return *this;
+ return static_cast<daughter &>(*this);
};
- inline iterator_internal & operator--() {
+ inline daughter & operator--() {
--ret;
- return *this;
+ return static_cast<daughter &>(*this);
};
- inline iterator_internal & operator+=(const UInt n) {
+ inline daughter & operator+=(const UInt n) {
ret += n;
- return *this;
+ return static_cast<daughter &>(*this);
}
- inline iterator_internal & operator-=(const UInt n) {
+ inline daughter & operator-=(const UInt n) {
ret -= n;
- return *this;
+ return static_cast<daughter &>(*this);
}
inline reference operator[](const UInt n) { return ret[n]; }
inline bool operator==(const iterator_internal & other) const {
return ret == other.ret;
}
inline bool operator!=(const iterator_internal & other) const {
return ret != other.ret;
}
inline bool operator<(const iterator_internal & other) const {
return ret < other.ret;
}
inline bool operator<=(const iterator_internal & other) const {
return ret <= other.ret;
}
inline bool operator>(const iterator_internal & other) const {
return ret > other.ret;
}
inline bool operator>=(const iterator_internal & other) const {
return ret >= other.ret;
}
- inline iterator_internal operator-(difference_type n) {
- return iterator_internal(ret - n);
- }
- inline iterator_internal operator+(difference_type n) {
- return iterator_internal(ret + n);
- }
+ inline daughter operator-(difference_type n) { return daughter(ret - n); }
+ inline daughter operator+(difference_type n) { return daughter(ret + n); }
inline difference_type operator-(const iterator_internal & b) {
return ret - b.ret;
}
inline pointer data() const { return ret; }
- inline difference_type offset() const { return _offset; }
protected:
- difference_type _offset;
- pointer ret;
- pointer initial;
+ pointer ret{nullptr};
+ pointer initial{nullptr};
};
/* -------------------------------------------------------------------------- */
/* Begin/End functions implementation */
/* -------------------------------------------------------------------------- */
-namespace {
- template <std::size_t N> struct extract_last {
- template <typename F, typename... T, typename... Arg>
- static decltype(auto) extract(F && func, std::tuple<T...> && t,
- Arg... args) {
- return extract_last<N - 1>::extract(
- std::forward<F>(func), std::forward<decltype(t)>(t),
- std::get<sizeof...(T) - N>(std::forward<decltype(t)>(t)), args...);
- }
- };
+namespace detail {
+ template <class Tuple, size_t... Is>
+ constexpr auto take_front_impl(Tuple && t, std::index_sequence<Is...>) {
+ return std::make_tuple(std::get<Is>(std::forward<Tuple>(t))...);
+ }
- template <> struct extract_last<1> {
- template <typename F, typename... T, typename... Arg>
- static decltype(auto) extract(F && func, std::tuple<T...> && /*unused*/,
- Arg... args) {
- return std::forward<F>(func)(args...);
- }
- };
+ template <size_t N, class Tuple> constexpr auto take_front(Tuple && t) {
+ return take_front_impl(std::forward<Tuple>(t),
+ std::make_index_sequence<N>{});
+ }
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wnarrowing"
- template <std::size_t N> struct InstantiationHelper {
- template <typename... Ns> static constexpr std::size_t product(Ns... ns) {
- std::size_t p = 1;
- for (auto n : std::array<std::size_t, sizeof...(Ns)>{ns...})
- p *= n;
- return p;
- }
+ template <typename... V>
+ constexpr auto product_all(V &&... v) ->
+ typename std::common_type<V...>::type {
+ typename std::common_type<V...>::type result = 1;
+ (void)std::initializer_list<int>{(result *= v, 0)...};
+ return result;
+ }
- template <typename... Ns> static std::string to_string(Ns... ns) {
- std::stringstream sstr;
- bool first = true;
- sstr << "(";
- for (auto n : std::array<std::size_t, sizeof...(Ns)>{ns...}) {
- if (!first) {
- sstr << ", ";
- }
- sstr << n;
- first = false;
- }
- sstr << ")";
- return sstr.str();
- }
-#pragma GCC diagnostic pop
+ template <typename... T> std::string to_string_all(T &&... t) {
+ if (sizeof...(T) == 0)
+ return "";
+
+ std::stringstream ss;
+ bool noComma = true;
+ ss << "(";
+ (void)std::initializer_list<bool>{
+ (ss << (noComma ? "" : ", ") << t, noComma = false)...};
+ ss << ")";
+ return ss.str();
+ }
+
+ template <std::size_t N> struct InstantiationHelper {
template <typename type, typename T, typename... Ns>
static auto instantiate(T && data, Ns... ns) {
- return new type(data, ns...);
+ return std::make_unique<type>(data, ns...);
}
};
template <> struct InstantiationHelper<0> {
template <typename type, typename T> static auto instantiate(T && data) {
return data;
}
-
- static constexpr std::size_t product() { return 1; }
- static std::string to_string() { return ""; }
};
template <typename Arr, typename T, typename... Ns>
- decltype(auto) get_iterator(Arr && array, T * data, Ns... ns) {
- using type = IteratorHelper_t<sizeof...(Ns) -1, T>;
+ decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) {
+ using type = IteratorHelper_t<sizeof...(Ns) - 1, T>;
using array_type = std::decay_t<Arr>;
using iterator =
std::conditional_t<std::is_const<Arr>::value,
typename array_type::template const_iterator<type>,
typename array_type::template iterator<type>>;
- AKANTU_DEBUG_ASSERT(
- array.getNbComponent() * array.size() ==
- InstantiationHelper<sizeof...(Ns)>::product(ns...),
- "The iterator is not compatible with the type "
- << debug::demangle(typeid(type).name())
- << InstantiationHelper<sizeof...(Ns)>::to_string(ns...));
+ if (array.getNbComponent() * array.size() !=
+ product_all(std::forward<Ns>(ns)...)) {
+ AKANTU_CUSTOM_EXCEPTION_INFO(
+ debug::ArrayException(),
+ "The iterator on Array "
+ << to_string_all(array.size(), array.getNbComponent())
+ << "is not compatible with the type "
+ << debug::demangle(typeid(type).name()) << to_string_all(ns...));
+ }
- auto && wrapped = extract_last<sizeof...(Ns)>::extract(
+ auto && wrapped = aka::apply(
[&](auto... n) {
return InstantiationHelper<sizeof...(n)>::template instantiate<type>(
data, n...);
},
- std::make_tuple(ns...));
+ take_front<sizeof...(Ns) - 1>(std::make_tuple(ns...)));
- return iterator{wrapped};
+ return iterator(std::move(wrapped));
}
-} // namespace
+} // namespace detail
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::begin(Ns... ns) {
- return get_iterator(*this, values, ns..., size_);
+inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) {
+ return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::end(Ns... ns) {
- return get_iterator(*this, values + nb_component * size_, ns..., size_);
+inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) {
+ return detail::get_iterator(*this, values + nb_component * size_,
+ std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::begin(Ns... ns) const {
- return get_iterator(*this, values, ns..., size_);
+inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) const {
+ return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::end(Ns... ns) const {
- return get_iterator(*this, values + nb_component * size_, ns..., size_);
+inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) const {
+ return detail::get_iterator(*this, values + nb_component * size_,
+ std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns... ns) {
- return get_iterator(*this, values, ns...);
+inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) {
+ return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns... ns) {
- return get_iterator(
- *this, values + InstantiationHelper<sizeof...(Ns)>::product(ns...),
- ns...);
+inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) {
+ return detail::get_iterator(
+ *this, values + detail::product_all(std::forward<Ns>(ns)...),
+ std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns... ns) const {
- return get_iterator(*this, values, ns...);
+inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) const {
+ return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
-inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns... ns) const {
- return get_iterator(
- *this, values + InstantiationHelper<sizeof...(Ns)>::product(ns...),
- ns...);
+inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) const {
+ return detail::get_iterator(
+ *this, values + detail::product_all(std::forward<Ns>(ns)...),
+ std::forward<Ns>(ns)...);
+}
+
+/* -------------------------------------------------------------------------- */
+/* Views */
+/* -------------------------------------------------------------------------- */
+namespace detail {
+ template <typename Array, typename... Ns> class ArrayView {
+ using tuple = std::tuple<Ns...>;
+
+ public:
+ ArrayView(Array && array, Ns &&... ns)
+ : array(std::forward<Array>(array)),
+ sizes(std::forward<Ns>(ns)...){};
+
+ decltype(auto) begin() {
+ return aka::apply(
+ [&](auto &&... ns) { return array.begin_reinterpret(ns...); }, sizes);
+ }
+
+ decltype(auto) end() {
+ return aka::apply(
+ [&](auto &&... ns) { return array.end_reinterpret(ns...); }, sizes);
+ }
+
+ decltype(auto) size() {
+ return std::get<std::tuple_size<tuple>::value - 1>(sizes);
+ }
+
+ private:
+ Array array;
+ tuple sizes;
+ };
+} // namespace detail
+
+/* -------------------------------------------------------------------------- */
+template <typename Array, typename... Ns>
+decltype(auto) make_view(Array && array, Ns &&... ns) {
+ auto size = std::forward<decltype(array)>(array).size() *
+ std::forward<decltype(array)>(array).getNbComponent() /
+ detail::product_all(ns...);
+
+ return detail::ArrayView<Array, Ns..., decltype(size)>(
+ std::forward<Array>(array), std::forward<Ns>(ns)...,
+ std::forward<decltype(size)>(size));
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
-class Array<T, is_scal>::const_iterator : public iterator_internal<const R, R> {
+class Array<T, is_scal>::const_iterator
+ : public iterator_internal<const R, Array<T, is_scal>::const_iterator<R>,
+ R> {
public:
- typedef iterator_internal<const R, R> parent;
+ typedef iterator_internal<const R, const_iterator, R> parent;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
const_iterator() : parent(){};
- const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}
- const_iterator(pointer warped) : parent(warped) {}
- const_iterator(const parent & it) : parent(it) {}
- // const_iterator(const const_iterator<R> & it) : parent(it) {}
+ // const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}
+ // const_iterator(pointer warped) : parent(warped) {}
+ // const_iterator(const parent & it) : parent(it) {}
- inline const_iterator operator+(difference_type n) {
- return parent::operator+(n);
- }
- inline const_iterator operator-(difference_type n) {
- return parent::operator-(n);
- }
- inline difference_type operator-(const const_iterator & b) {
- return parent::operator-(b);
- }
+ const_iterator(const const_iterator & it) = default;
+ const_iterator(const_iterator && it) = default;
- inline const_iterator & operator++() {
- parent::operator++();
- return *this;
- };
- inline const_iterator & operator--() {
- parent::operator--();
- return *this;
- };
- inline const_iterator & operator+=(const UInt n) {
- parent::operator+=(n);
- return *this;
- }
+ template <typename P, typename = std::enable_if_t<not is_tensor<P>{}>>
+ const_iterator(P * data) : parent(data) {}
+
+ template <typename UP_P, typename = std::enable_if_t<
+ is_tensor<typename UP_P::element_type>::value>>
+ const_iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
+
+ const_iterator & operator=(const const_iterator & it) = default;
};
-// #endif
-// #if defined(AKANTU_CORE_CXX11)
-// template<class R> using iterator = iterator_internal<R>;
-// #else
-template <class T, class R, bool issame = std::is_same<T, R>::value>
+template <class T, class R, bool is_tensor_ = is_tensor<R>{}>
struct ConstConverterIteratorHelper {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
+
static inline const_iterator convert(const iterator & it) {
- return const_iterator(new R(*it, false));
+ return const_iterator(std::unique_ptr<R>(new R(*it, false)));
}
};
-template <class T, class R> struct ConstConverterIteratorHelper<T, R, true> {
+template <class T, class R> struct ConstConverterIteratorHelper<T, R, false> {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
static inline const_iterator convert(const iterator & it) {
- return const_iterator(it.data(), it.offset());
+ return const_iterator(it.data());
}
};
template <class T, bool is_scal>
template <typename R>
-class Array<T, is_scal>::iterator : public iterator_internal<R> {
+class Array<T, is_scal>::iterator
+ : public iterator_internal<R, Array<T, is_scal>::iterator<R>> {
public:
- using parent = iterator_internal<R>;
+ using parent = iterator_internal<R, iterator>;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
iterator() : parent(){};
- iterator(pointer_type data, UInt offset) : parent(data, offset){};
- iterator(pointer warped) : parent(warped) {}
- iterator(const parent & it) : parent(it) {}
- // iterator(const iterator<R> & it) : parent(it) {}
+ iterator(const iterator & it) = default;
+ iterator(iterator && it) = default;
- operator const_iterator<R>() {
- return ConstConverterIteratorHelper<T, R>::convert(*this);
- }
+ template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
+ iterator(P * data) : parent(data) {}
- inline iterator operator+(difference_type n) {
- return parent::operator+(n);
- ;
- }
- inline iterator operator-(difference_type n) {
- return parent::operator-(n);
- ;
- }
- inline difference_type operator-(const iterator & b) {
- return parent::operator-(b);
- }
+ template <typename UP_P, typename = std::enable_if_t<
+ is_tensor<typename UP_P::element_type>::value>>
+ iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
- inline iterator & operator++() {
- parent::operator++();
- return *this;
- };
- inline iterator & operator--() {
- parent::operator--();
- return *this;
- };
- inline iterator & operator+=(const UInt n) {
- parent::operator+=(n);
- return *this;
+ iterator & operator=(const iterator & it) = default;
+
+ operator const_iterator<R>() {
+ return ConstConverterIteratorHelper<T, R>::convert(*this);
}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
inline typename Array<T, is_scal>::template iterator<R>
Array<T, is_scal>::erase(const iterator<R> & it) {
T * curr = it.data();
UInt pos = (curr - values) / nb_component;
erase(pos);
iterator<R> rit = it;
return --rit;
}
} // namespace akantu
#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */
diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh
index b21921277..50c1bcc8a 100644
--- a/src/common/aka_common.hh
+++ b/src/common/aka_common.hh
@@ -1,415 +1,414 @@
/**
* @file aka_common.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Thu Jan 21 2016
*
* @brief common type descriptions for akantu
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 <http://www.gnu.org/licenses/>.
*
* @section DESCRIPTION
*
* All common things to be included in the projects files
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COMMON_HH__
#define __AKANTU_COMMON_HH__
/* -------------------------------------------------------------------------- */
-#include <list>
#include <limits>
+#include <list>
#include <type_traits>
/* -------------------------------------------------------------------------- */
#include "aka_compatibilty_with_cpp_standard.hh"
/* -------------------------------------------------------------------------- */
#define __BEGIN_AKANTU_DUMPER__ namespace dumper {
#define __END_AKANTU_DUMPER__ }
/* -------------------------------------------------------------------------- */
#if defined(WIN32)
#define __attribute__(x)
#endif
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
#include "aka_error.hh"
#include "aka_safe_enum.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Common types */
/* -------------------------------------------------------------------------- */
using ID = std::string;
#ifdef AKANTU_NDEBUG
static const Real REAL_INIT_VALUE = Real(0.);
#else
static const Real REAL_INIT_VALUE = std::numeric_limits<Real>::quiet_NaN();
#endif
/* -------------------------------------------------------------------------- */
/* Memory types */
/* -------------------------------------------------------------------------- */
using MemoryID = UInt;
using Surface = std::string;
typedef std::pair<Surface, Surface> SurfacePair;
using SurfacePairList = std::list<SurfacePair>;
/* -------------------------------------------------------------------------- */
extern const UInt _all_dimensions;
/* -------------------------------------------------------------------------- */
/* Mesh/FEM/Model types */
/* -------------------------------------------------------------------------- */
-} // akantu
+} // namespace akantu
#include "aka_element_classes_info.hh"
namespace akantu {
/// small help to use names for directions
enum SpacialDirection { _x = 0, _y = 1, _z = 2 };
/// enum MeshIOType type of mesh reader/writer
enum MeshIOType {
_miot_auto, ///< Auto guess of the reader to use based on the extension
_miot_gmsh, ///< Gmsh files
_miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has
/// structures elements
_miot_diana, ///< TNO Diana mesh format
_miot_abaqus ///< Abaqus mesh format
};
/// enum MeshEventHandlerPriority defines relative order of execution of events
enum EventHandlerPriority {
_ehp_highest = 0,
_ehp_mesh = 5,
_ehp_fe_engine = 9,
_ehp_synchronizer = 10,
_ehp_dof_manager = 20,
_ehp_model = 94,
_ehp_non_local_manager = 100,
_ehp_lowest = 100
};
-
/// enum AnalysisMethod type of solving method used to solve the equation of
/// motion
enum AnalysisMethod {
_static = 0,
_implicit_dynamic = 1,
_explicit_lumped_mass = 2,
_explicit_lumped_capacity = 2,
_explicit_consistent_mass = 3
};
/// enum DOFSupportType defines which kind of dof that can exists
enum DOFSupportType { _dst_nodal, _dst_generic };
/// Type of non linear resolution available in akantu
enum NonLinearSolverType {
_nls_linear, ///< No non linear convergence loop
_nls_newton_raphson, ///< Regular Newton-Raphson
_nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent
_nls_lumped, ///< Case of lumped mass or equivalent matrix
- _nls_auto ///< This will take a default value that make sense in case of
- /// model::getNewSolver
+ _nls_auto ///< This will take a default value that make sense in case of
+ /// model::getNewSolver
};
/// Define the node/dof type
-enum NodeType : Int {
- _nt_pure_gost = -3,
- _nt_master = -2,
- _nt_normal = -1
-};
+enum NodeType : Int { _nt_pure_gost = -3, _nt_master = -2, _nt_normal = -1 };
/// Type of time stepping solver
enum TimeStepSolverType {
_tsst_static, ///< Static solution
_tsst_dynamic, ///< Dynamic solver
_tsst_dynamic_lumped, ///< Dynamic solver with lumped mass
_tsst_not_defined, ///< For not defined cases
};
/// Type of integration scheme
enum IntegrationSchemeType {
- _ist_pseudo_time, ///< Pseudo Time
- _ist_forward_euler, ///< GeneralizedTrapezoidal(0)
- _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2)
- _ist_backward_euler, ///< GeneralizedTrapezoidal(1)
- _ist_central_difference, ///< NewmarkBeta(0, 1/2)
- _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2)
- _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2)
- _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2)
- _ist_newmark_beta, ///< generic NewmarkBeta with user defined
- /// alpha and beta
- _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user
- /// defined alpha
+ _ist_pseudo_time, ///< Pseudo Time
+ _ist_forward_euler, ///< GeneralizedTrapezoidal(0)
+ _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2)
+ _ist_backward_euler, ///< GeneralizedTrapezoidal(1)
+ _ist_central_difference, ///< NewmarkBeta(0, 1/2)
+ _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2)
+ _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2)
+ _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2)
+ _ist_newmark_beta, ///< generic NewmarkBeta with user defined
+ /// alpha and beta
+ _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user
+ /// defined alpha
};
/// enum SolveConvergenceCriteria different convergence criteria
enum SolveConvergenceCriteria {
_scc_residual, ///< Use residual to test the convergence
_scc_solution, ///< Use solution to test the convergence
_scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb
};
/// enum CohesiveMethod type of insertion of cohesive elements
enum CohesiveMethod { _intrinsic, _extrinsic };
/// @enum SparseMatrixType type of sparse matrix used
enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined };
/* -------------------------------------------------------------------------- */
/* Ghosts handling */
/* -------------------------------------------------------------------------- */
/// @enum CommunicatorType type of communication method to use
enum CommunicatorType { _communicator_mpi, _communicator_dummy };
/// @enum SynchronizationTag type of synchronizations
enum SynchronizationTag {
//--- Generic tags ---
_gst_whatever,
_gst_update,
_gst_size,
//--- SolidMechanicsModel tags ---
_gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass
_gst_smm_for_gradu, ///< synchronization of the
/// SolidMechanicsModel.displacement
_gst_smm_boundary, ///< synchronization of the boundary, forces, velocities
/// and displacement
_gst_smm_uv, ///< synchronization of the nodal velocities and displacement
_gst_smm_res, ///< synchronization of the nodal residual
_gst_smm_init_mat, ///< synchronization of the data to initialize materials
_gst_smm_stress, ///< synchronization of the stresses to compute the internal
/// forces
_gst_smmc_facets, ///< synchronization of facet data to setup facet synch
_gst_smmc_facets_conn, ///< synchronization of facet global connectivity
_gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet
/// synch
_gst_smmc_damage, ///< synchronization of damage
// --- GlobalIdsUpdater tags ---
_gst_giu_global_conn, ///< synchronization of global connectivities
// --- CohesiveElementInserter tags ---
_gst_ce_groups, ///< synchronization of cohesive element insertion depending
/// on facet groups
// --- GroupManager tags ---
_gst_gm_clusters, ///< synchronization of clusters
// --- HeatTransfer tags ---
_gst_htm_capacity, ///< synchronization of the nodal heat capacity
_gst_htm_temperature, ///< synchronization of the nodal temperature
_gst_htm_gradient_temperature, ///< synchronization of the element gradient
/// temperature
// --- LevelSet tags ---
_gst_htm_phi, ///< synchronization of the nodal level set value phi
_gst_htm_gradient_phi, ///< synchronization of the element gradient phi
//--- Material non local ---
_gst_mnl_for_average, ///< synchronization of data to average in non local
/// material
_gst_mnl_weight, ///< synchronization of data for the weight computations
// --- NeighborhoodSynchronization tags ---
_gst_nh_criterion,
// --- General tags ---
_gst_test, ///< Test tag
_gst_user_1, ///< tag for user simulations
_gst_user_2, ///< tag for user simulations
_gst_material_id, ///< synchronization of the material ids
_gst_for_dump, ///< everything that needs to be synch before dump
// --- Contact & Friction ---
_gst_cf_nodal, ///< synchronization of disp, velo, and current position
_gst_cf_incr, ///< synchronization of increment
// --- Solver tags ---
_gst_solver_solution ///< synchronization of the solution obained with the
/// PETSc solver
};
/// standard output stream operator for SynchronizationTag
inline std::ostream & operator<<(std::ostream & stream,
SynchronizationTag type);
/// @enum GhostType type of ghost
enum GhostType {
_not_ghost,
_ghost,
_casper // not used but a real cute ghost
};
/* -------------------------------------------------------------------------- */
struct GhostType_def {
using type = GhostType;
static const type _begin_ = _not_ghost;
static const type _end_ = _casper;
};
using ghost_type_t = safe_enum<GhostType_def>;
extern ghost_type_t ghost_types;
/// standard output stream operator for GhostType
inline std::ostream & operator<<(std::ostream & stream, GhostType type);
/* -------------------------------------------------------------------------- */
/* Global defines */
/* -------------------------------------------------------------------------- */
#define AKANTU_MIN_ALLOCATION 2000
#define AKANTU_INDENT " "
#define AKANTU_INCLUDE_INLINE_IMPL
/* -------------------------------------------------------------------------- */
-template<typename T>
-using is_scalar = std::is_arithmetic<T>;
+/* Type traits */
+/* -------------------------------------------------------------------------- */
+struct TensorTrait {};
+/* -------------------------------------------------------------------------- */
+template <typename T> using is_tensor = std::is_base_of<TensorTrait, T>;
+/* -------------------------------------------------------------------------- */
+template <typename T> using is_scalar = std::is_arithmetic<T>;
+/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#define AKANTU_SET_MACRO(name, variable, type) \
inline void set##name(type variable) { this->variable = variable; }
#define AKANTU_GET_MACRO(name, variable, type) \
inline type get##name() const { return variable; }
#define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \
inline type get##name() { return variable; }
#define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \
inline con Array<type> & get##name( \
const support & el_type, const GhostType & ghost_type = _not_ghost) \
con { \
return variable(el_type, ghost_type); \
}
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, )
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const)
#define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, )
#define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const)
/* -------------------------------------------------------------------------- */
/// initialize the static part of akantu
void initialize(int & argc, char **& argv);
/// initialize the static part of akantu and read the global input_file
void initialize(const std::string & input_file, int & argc, char **& argv);
/* -------------------------------------------------------------------------- */
/// finilize correctly akantu and clean the memory
void finalize();
/* -------------------------------------------------------------------------- */
/// Read an new input file
void readInputFile(const std::string & input_file);
/* -------------------------------------------------------------------------- */
/*
* For intel compiler annoying remark
*/
// #if defined(__INTEL_COMPILER)
// /// remark #981: operands are evaluated in unspecified order
// #pragma warning(disable : 981)
// /// remark #383: value copied to temporary, reference to temporary used
// #pragma warning(disable : 383)
// #endif // defined(__INTEL_COMPILER)
/* -------------------------------------------------------------------------- */
/* string manipulation */
/* -------------------------------------------------------------------------- */
inline std::string to_lower(const std::string & str);
/* -------------------------------------------------------------------------- */
inline std::string trim(const std::string & to_trim);
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/// give a string representation of the a human readable size in bit
template <typename T> std::string printMemorySize(UInt size);
/* -------------------------------------------------------------------------- */
-} // akantu
+} // namespace akantu
#include "aka_fwd.hh"
namespace akantu {
/// get access to the internal argument parser
cppargparse::ArgumentParser & getStaticArgumentParser();
/// get access to the internal input file parser
Parser & getStaticParser();
/// get access to the user part of the internal input file parser
const ParserSection & getUserParser();
-} // akantu
+} // namespace akantu
#include "aka_common_inline_impl.cc"
/* -------------------------------------------------------------------------- */
#if AKANTU_INTEGER_SIZE == 4
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9
#elif AKANTU_INTEGER_SIZE == 8
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL
#endif
namespace std {
/**
* Hashing function for pairs based on hash_combine from boost The magic number
* is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f]
* @f[\frac{2^32}{\phi} = 0x9e3779b9@f]
* http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine
* http://burtleburtle.net/bob/hash/doobs.html
*/
-template <typename a, typename b> struct hash<std::pair<a, b> > {
-public:
- hash() : ah(), bh() {}
+template <typename a, typename b> struct hash<std::pair<a, b>> {
+ hash() = default;
size_t operator()(const std::pair<a, b> & p) const {
size_t seed = ah(p.first);
return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) +
(seed >> 2);
}
private:
- const hash<a> ah;
- const hash<b> bh;
+ const hash<a> ah{};
+ const hash<b> bh{};
};
-} //std
-
+} // namespace std
#endif /* __AKANTU_COMMON_HH__ */
diff --git a/src/common/aka_compatibilty_with_cpp_standard.hh b/src/common/aka_compatibilty_with_cpp_standard.hh
index c3b223671..154f339f1 100644
--- a/src/common/aka_compatibilty_with_cpp_standard.hh
+++ b/src/common/aka_compatibilty_with_cpp_standard.hh
@@ -1,63 +1,132 @@
/**
* @file aka_compatibilty_with_cpp_standard.hh
*
* @author Nicolas Richart
*
* @date creation Wed Oct 25 2017
*
* @brief The content of this file is taken from the possible implementations on
* http://en.cppreference.com
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <type_traits>
+#include <utility>
+#include <tuple>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_COMPATIBILTY_WITH_CPP_STANDARD_HH__
#define __AKANTU_AKA_COMPATIBILTY_WITH_CPP_STANDARD_HH__
-namespace std {
+namespace aka {
// Part taken from C++14
#if __cplusplus < 201402L
template <bool B, class T = void>
using enable_if_t = typename enable_if<B, T>::type;
+#else
+template <bool B, class T = void>
+using enable_if_t = std::enable_if_t<B, T>;
#endif
// Part taken from C++17
#if __cplusplus < 201703L
// bool_constant
-template <bool B> using bool_constant = integral_constant<bool, B>;
+template <bool B> using bool_constant = std::integral_constant<bool, B>;
namespace {
template <bool B> constexpr bool bool_constant_v = bool_constant<B>::value;
}
// conjunction
template <class...> struct conjunction : std::true_type {};
template <class B1> struct conjunction<B1> : B1 {};
template <class B1, class... Bn>
struct conjunction<B1, Bn...>
: std::conditional_t<bool(B1::value), conjunction<Bn...>, B1> {};
-#endif
+namespace detail {
+ // template <class T>
+ // struct is_reference_wrapper : std::false_type {};
+ // template <class U>
+ // struct is_reference_wrapper<std::reference_wrapper<U>> : std::true_type {};
+ // template <class T>
+ // constexpr bool is_reference_wrapper_v = is_reference_wrapper<T>::value;
+
+ template <class T, class Type, class T1, class... Args>
+ decltype(auto) INVOKE(Type T::*f, T1 && t1, Args &&... args) {
+ static_assert(std::is_member_function_pointer<decltype(f)>{} and
+ std::is_base_of<T, std::decay_t<T1>>{},
+ "Does not know what to do with this types");
+ return (std::forward<T1>(t1).*f)(std::forward<Args>(args)...);
+ }
+
+ // template <class T, class Type, class T1, class... Args>
+ // decltype(auto) INVOKE(Type T::*f, T1 && t1, Args &&... args) {
+ // if constexpr (std::is_member_function_pointer_v<decltype(f)>) {
+ // if constexpr (std::is_base_of_v<T, std::decay_t<T1>>)
+ // return (std::forward<T1>(t1).*f)(std::forward<Args>(args)...);
+ // else if constexpr (is_reference_wrapper_v<std::decay_t<T1>>)
+ // return (t1.get().*f)(std::forward<Args>(args)...);
+ // else
+ // return ((*std::forward<T1>(t1)).*f)(std::forward<Args>(args)...);
+ // } else {
+ // static_assert(std::is_member_object_pointer_v<decltype(f)>);
+ // static_assert(sizeof...(args) == 0);
+ // if constexpr (std::is_base_of_v<T, std::decay_t<T1>>)
+ // return std::forward<T1>(t1).*f;
+ // else if constexpr (is_reference_wrapper_v<std::decay_t<T1>>)
+ // return t1.get().*f;
+ // else
+ // return (*std::forward<T1>(t1)).*f;
+ // }
+ // }
+
+ template <class F, class... Args>
+ decltype(auto) INVOKE(F && f, Args &&... args) {
+ return std::forward<F>(f)(std::forward<Args>(args)...);
+ }
+} // namespace detail
+
+template <class F, class... Args>
+decltype(auto) invoke(F && f, Args &&... args) {
+ return detail::INVOKE(std::forward<F>(f), std::forward<Args>(args)...);
}
+namespace detail {
+ template <class F, class Tuple, std::size_t... Is>
+ constexpr decltype(auto) apply_impl(F && f, Tuple && t,
+ std::index_sequence<Is...>) {
+ return invoke(std::forward<F>(f),
+ std::get<Is>(std::forward<Tuple>(t))...);
+ }
+} // namespace detail
+
+template <class F, class Tuple>
+constexpr decltype(auto) apply(F && f, Tuple && t) {
+ return detail::apply_impl(
+ std::forward<F>(f), std::forward<Tuple>(t),
+ std::make_index_sequence<std::tuple_size<std::decay_t<Tuple>>::value>{});
+}
+
+#endif
+} // namespace aka
+
#endif /* __AKANTU_AKA_COMPATIBILTY_WITH_CPP_STANDARD_HH__ */
diff --git a/src/common/aka_iterators.hh b/src/common/aka_iterators.hh
index e5a41ca92..f76cf6214 100644
--- a/src/common/aka_iterators.hh
+++ b/src/common/aka_iterators.hh
@@ -1,326 +1,341 @@
/**
* @file aka_iterators.hh
*
* @author Nicolas Richart
*
* @date creation Wed Jul 19 2017
*
* @brief iterator interfaces
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_error.hh"
/* -------------------------------------------------------------------------- */
#include <cassert>
#include <iostream>
#include <tuple>
#include <utility>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_ITERATORS_HH__
#define __AKANTU_AKA_ITERATORS_HH__
namespace akantu {
namespace tuple {
/* ------------------------------------------------------------------------ */
namespace details {
template <size_t N> struct Foreach {
template <class F, class Tuple>
static inline decltype(auto) transform_forward(F && func,
Tuple && tuple) {
return std::tuple_cat(
Foreach<N - 1>::transform_forward(std::forward<F>(func),
std::forward<Tuple>(tuple)),
std::forward_as_tuple(std::forward<F>(func)(
std::get<N - 1>(std::forward<Tuple>(tuple)))));
}
template <class F, class Tuple>
static inline decltype(auto) transform(F && func, Tuple && tuple) {
return std::tuple_cat(
Foreach<N - 1>::transform(std::forward<F>(func),
std::forward<Tuple>(tuple)),
std::make_tuple(std::forward<F>(func)(
std::get<N - 1>(std::forward<Tuple>(tuple)))));
}
template <class F, class Tuple>
static inline void foreach (F && func, Tuple && tuple) {
Foreach<N - 1>::foreach (std::forward<F>(func),
std::forward<Tuple>(tuple));
std::forward<F>(func)(std::get<N - 1>(std::forward<Tuple>(tuple)));
}
template <class Tuple> static inline bool equal(Tuple && a, Tuple && b) {
if (not(std::get<N - 1>(std::forward<Tuple>(a)) ==
std::get<N - 1>(std::forward<Tuple>(b))))
return false;
return Foreach<N - 1>::equal(std::forward<Tuple>(a),
std::forward<Tuple>(b));
}
};
/* ------------------------------------------------------------------------
*/
template <> struct Foreach<1> {
template <class F, class Tuple>
static inline decltype(auto) transform_forward(F && func,
Tuple && tuple) {
return std::forward_as_tuple(
std::forward<F>(func)(std::get<0>(std::forward<Tuple>(tuple))));
}
template <class F, class Tuple>
static inline decltype(auto) transform(F && func, Tuple && tuple) {
return std::make_tuple(
std::forward<F>(func)(std::get<0>(std::forward<Tuple>(tuple))));
}
template <class F, class Tuple>
static inline void foreach (F && func, Tuple && tuple) {
std::forward<F>(func)(std::get<0>(std::forward<Tuple>(tuple)));
}
template <class Tuple> static inline bool equal(Tuple && a, Tuple && b) {
return std::get<0>(std::forward<Tuple>(a)) ==
std::get<0>(std::forward<Tuple>(b));
}
};
} // namespace details
/* ------------------------------------------------------------------------ */
template <class Tuple> bool are_equal(Tuple && a, Tuple && b) {
return details::Foreach<std::tuple_size<std::decay_t<Tuple>>::value>::equal(
std::forward<Tuple>(a), std::forward<Tuple>(b));
}
template <class F, class Tuple> void foreach (F && func, Tuple && tuple) {
details::Foreach<std::tuple_size<std::decay_t<Tuple>>::value>::foreach (
std::forward<F>(func), std::forward<Tuple>(tuple));
}
template <class F, class Tuple>
decltype(auto) transform_forward(F && func, Tuple && tuple) {
return details::Foreach<std::tuple_size<std::decay_t<Tuple>>::value>::
transform_forward(std::forward<F>(func), std::forward<Tuple>(tuple));
}
template <class F, class Tuple>
decltype(auto) transform(F && func, Tuple && tuple) {
return details::Foreach<std::tuple_size<std::decay_t<Tuple>>::value>::
transform(std::forward<F>(func), std::forward<Tuple>(tuple));
}
} // namespace tuple
namespace iterators {
namespace details {
struct dereference_iterator {
template <class Iter> decltype(auto) operator()(Iter & it) const {
return std::forward<decltype(*it)>(*it);
}
};
struct increment_iterator {
template <class Iter> void operator()(Iter & it) const { ++it; }
};
struct begin_container {
template <class Container>
decltype(auto) operator()(Container && cont) const {
return std::forward<Container>(cont).begin();
}
};
struct end_container {
template <class Container>
decltype(auto) operator()(Container && cont) const {
return std::forward<Container>(cont).end();
}
};
} // namespace details
/* ------------------------------------------------------------------------ */
template <class... Iterators> class ZipIterator {
private:
using tuple_t = std::tuple<Iterators...>;
public:
explicit ZipIterator(tuple_t iterators) : iterators(std::move(iterators)) {}
decltype(auto) operator*() {
return tuple::transform_forward(details::dereference_iterator(),
iterators);
}
ZipIterator & operator++() {
tuple::foreach (details::increment_iterator(), iterators);
return *this;
}
bool operator==(const ZipIterator & other) const {
return tuple::are_equal(iterators, other.iterators);
}
bool operator!=(const ZipIterator & other) const {
return not operator==(other);
}
private:
tuple_t iterators;
};
} // namespace iterators
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <class... Iterators>
decltype(auto) zip_iterator(std::tuple<Iterators...> && iterators_tuple) {
auto zip = iterators::ZipIterator<Iterators...>(
std::forward<decltype(iterators_tuple)>(iterators_tuple));
return zip;
}
/* -------------------------------------------------------------------------- */
namespace containers {
template <class... Containers> class ZipContainer {
using containers_t = std::tuple<Containers...>;
public:
explicit ZipContainer(Containers &&... containers)
: containers(std::forward<Containers>(containers)...) {}
decltype(auto) begin() const {
return zip_iterator(
tuple::transform(iterators::details::begin_container(),
std::forward<containers_t>(containers)));
}
decltype(auto) end() const {
return zip_iterator(
tuple::transform(iterators::details::end_container(),
std::forward<containers_t>(containers)));
}
decltype(auto) begin() {
return zip_iterator(
tuple::transform(iterators::details::begin_container(),
std::forward<containers_t>(containers)));
}
decltype(auto) end() {
return zip_iterator(
tuple::transform(iterators::details::end_container(),
std::forward<containers_t>(containers)));
}
private:
containers_t containers;
};
} // namespace containers
/* -------------------------------------------------------------------------- */
template <class... Containers> decltype(auto) zip(Containers &&... conts) {
return containers::ZipContainer<Containers...>(
std::forward<Containers>(conts)...);
}
/* -------------------------------------------------------------------------- */
/* Arange */
/* -------------------------------------------------------------------------- */
namespace iterators {
template <class T> class ArangeIterator {
public:
using value_type = T;
using pointer = T *;
using reference = T &;
using iterator_category = std::input_iterator_tag;
constexpr ArangeIterator(T value, T step) : value(value), step(step) {}
constexpr ArangeIterator(const ArangeIterator &) = default;
constexpr ArangeIterator & operator++() {
value += step;
return *this;
}
constexpr const T & operator*() const { return value; }
constexpr bool operator==(const ArangeIterator & other) const {
return (value == other.value) and (step == other.step);
}
constexpr bool operator!=(const ArangeIterator & other) const {
return not operator==(other);
}
private:
T value{0};
const T step{1};
};
} // namespace iterators
namespace containers {
template <class T> class ArangeContainer {
public:
using iterator = iterators::ArangeIterator<T>;
constexpr ArangeContainer(T start, T stop, T step = 1)
: start(start), stop((stop - start) % step == 0
? stop
: start + (1 + (stop - start) / step) * step),
step(step) {}
explicit constexpr ArangeContainer(T stop) : ArangeContainer(0, stop, 1) {}
constexpr T operator[](size_t i) {
T val = start + i * step;
assert(val < stop && "i is out of range");
return val;
}
- constexpr size_t size() { return (stop - start) / step; }
+ constexpr T size() { return (stop - start) / step; }
constexpr iterator begin() { return iterator(start, step); }
constexpr iterator end() { return iterator(stop, step); }
private:
const T start{0}, stop{0}, step{1};
};
} // namespace containers
-template <class T> inline decltype(auto) arange(T stop) {
+template <class T,
+ typename = std::enable_if_t<std::is_integral<std::decay_t<T>>::value>>
+inline decltype(auto) arange(const T & stop) {
return containers::ArangeContainer<T>(stop);
}
-template <class T>
-inline constexpr decltype(auto) arange(T start, T stop, T step = 1) {
- return containers::ArangeContainer<T>(start, stop, step);
+template <class T1, class T2,
+ typename = std::enable_if_t<
+ std::is_integral<std::common_type_t<T1, T2>>::value>>
+inline constexpr decltype(auto) arange(const T1 & start, const T2 & stop) {
+ return containers::ArangeContainer<std::common_type_t<T1, T2>>(start, stop);
}
+template <class T1, class T2, class T3,
+ typename = std::enable_if_t<
+ std::is_integral<std::common_type_t<T1, T2, T3>>::value>>
+inline constexpr decltype(auto) arange(const T1 & start, const T2 & stop,
+ const T3 & step) {
+ return containers::ArangeContainer<std::common_type_t<T1, T2, T3>>(
+ start, stop, step);
+}
+
+/* -------------------------------------------------------------------------- */
+
template <class Container>
inline constexpr decltype(auto) enumerate(Container && container,
size_t start_ = 0) {
auto stop = std::forward<Container>(container).size();
decltype(stop) start = start_;
return zip(arange(start, stop), std::forward<Container>(container));
}
} // namespace akantu
#endif /* __AKANTU_AKA_ITERATORS_HH__ */
diff --git a/src/common/aka_safe_enum.hh b/src/common/aka_safe_enum.hh
index 20af59899..b129b7fbd 100644
--- a/src/common/aka_safe_enum.hh
+++ b/src/common/aka_safe_enum.hh
@@ -1,83 +1,86 @@
/**
* @file aka_safe_enum.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 21 2013
* @date last modification: Sun Oct 19 2014
*
* @brief Safe enums type (see More C++ Idioms/Type Safe Enum on Wikibooks
* http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Type_Safe_Enum)
*
* @section LICENSE
*
* Copyright (©) 2014, 2015 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_SAFE_ENUM_HH__
#define __AKANTU_AKA_SAFE_ENUM_HH__
namespace akantu {
/// Safe enumerated type
template<typename def, typename inner = typename def::type>
class safe_enum : public def {
using type = typename def::type;
public:
explicit safe_enum(type v = def::_end_) : val(v) {}
+ safe_enum(safe_enum && other) = default;
+ safe_enum& operator=(safe_enum && other) = default;
+
inner underlying() const { return val; }
bool operator == (const safe_enum & s) const { return this->val == s.val; }
bool operator != (const safe_enum & s) const { return this->val != s.val; }
bool operator < (const safe_enum & s) const { return this->val < s.val; }
bool operator <= (const safe_enum & s) const { return this->val <= s.val; }
bool operator > (const safe_enum & s) const { return this->val > s.val; }
bool operator >= (const safe_enum & s) const { return this->val >= s.val; }
operator inner() { return val; };
public:
// Works only if enumerations are contiguous.
class iterator {
public:
explicit iterator(type v) : it(v) { }
iterator & operator++() { ++it; return *this; }
safe_enum operator*() { return safe_enum(static_cast<type>(it)); }
bool operator!=(iterator const & it) { return it.it != this->it; }
private:
int it;
};
static iterator begin() {
return iterator(def::_begin_);
}
static iterator end() {
return iterator(def::_end_);
}
protected:
inner val;
};
} // akantu
#endif /* __AKANTU_AKA_SAFE_ENUM_HH__ */
diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh
index d8da67f62..bd65ef234 100644
--- a/src/common/aka_types.hh
+++ b/src/common/aka_types.hh
@@ -1,1255 +1,1263 @@
/**
* @file aka_types.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 17 2011
* @date last modification: Fri Jan 22 2016
*
* @brief description of the "simple" types
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_error.hh"
#include "aka_fwd.hh"
#include "aka_math.hh"
/* -------------------------------------------------------------------------- */
#include <initializer_list>
#include <iomanip>
#include <type_traits>
/* -------------------------------------------------------------------------- */
#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 <UInt dim> 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 <typename T, UInt ndim, class RetType> class TensorStorage;
/* -------------------------------------------------------------------------- */
/* Proxy classes */
/* -------------------------------------------------------------------------- */
namespace tensors {
-template <class A, class B> struct is_copyable {
- enum : bool { value = false };
-};
+ template <class A, class B> struct is_copyable {
+ enum : bool { value = false };
+ };
-template <class A> struct is_copyable<A, A> {
- enum : bool { value = true };
-};
+ template <class A> struct is_copyable<A, A> {
+ enum : bool { value = true };
+ };
-template <class A> struct is_copyable<A, typename A::RetType> {
- enum : bool { value = true };
-};
+ template <class A> struct is_copyable<A, typename A::RetType> {
+ enum : bool { value = true };
+ };
-template <class A> struct is_copyable<A, typename A::RetType::proxy> {
- enum : bool { value = true };
-};
+ template <class A> struct is_copyable<A, typename A::RetType::proxy> {
+ 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<T>::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 <typename T, UInt ndim, class _RetType> class TensorProxy {
protected:
using RetTypeProxy = typename _RetType::proxy;
TensorProxy(T * data, UInt m, UInt n, UInt p) {
DimHelper<ndim>::setDims(m, n, p, this->n);
this->values = data;
}
template <class Other, typename = std::enable_if_t<
tensors::is_copyable<TensorProxy, Other>::value>>
explicit TensorProxy(const Other & other) {
this->values = other.storage();
for (UInt i = 0; i < ndim; ++i)
this->n[i] = other.size(i);
}
public:
using RetType = _RetType;
- operator RetType() { return RetType(static_cast<RetTypeProxy &>(*this)); }
+ //operator RetType() { return RetType(static_cast<RetTypeProxy &>(*this)); }
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; }
template <class Other, typename = std::enable_if_t<
tensors::is_copyable<TensorProxy, Other>::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 <class Other, typename = std::enable_if_t<
// tensors::is_copyable<TensorProxy, Other>::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 <typename O> inline RetTypeProxy & operator*=(const O & o) {
RetType(*this) *= o;
return static_cast<RetTypeProxy &>(*this);
}
template <typename O> inline RetTypeProxy & operator/=(const O & o) {
RetType(*this) /= o;
return static_cast<RetTypeProxy &>(*this);
}
protected:
T * values;
UInt n[ndim];
};
/* -------------------------------------------------------------------------- */
template <typename T> class VectorProxy : public TensorProxy<T, 1, Vector<T>> {
using parent = TensorProxy<T, 1, Vector<T>>;
using type = Vector<T>;
public:
VectorProxy(T * data, UInt n) : parent(data, n, 0, 0) {}
template <class Other> explicit VectorProxy(Other & src) : parent(src) {}
/* ---------------------------------------------------------------------- */
template <class Other>
inline VectorProxy<T> & operator=(const Other & other) {
parent::operator=(other);
return *this;
}
// inline VectorProxy<T> & 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 <typename T> class MatrixProxy : public TensorProxy<T, 2, Matrix<T>> {
using parent = TensorProxy<T, 2, Matrix<T>>;
using type = Matrix<T>;
public:
MatrixProxy(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
template <class Other> explicit MatrixProxy(Other & src) : parent(src) {}
/* ---------------------------------------------------------------------- */
template <class Other>
inline MatrixProxy<T> & operator=(const Other & other) {
parent::operator=(other);
return *this;
}
};
template <typename T>
class Tensor3Proxy : public TensorProxy<T, 3, Tensor3<T>> {
typedef TensorProxy<T, 3, Tensor3<T>> parent;
using type = Tensor3<T>;
public:
Tensor3Proxy(T * data, UInt m, UInt n, UInt k) : parent(data, m, n, k) {}
Tensor3Proxy(const Tensor3Proxy & src) : parent(src) {}
Tensor3Proxy(const Tensor3<T> & src) : parent(src) {}
/* ---------------------------------------------------------------------- */
template <class Other>
inline Tensor3Proxy<T> & operator=(const Other & other) {
parent::operator=(other);
return *this;
}
};
/* -------------------------------------------------------------------------- */
/* Tensor base class */
/* -------------------------------------------------------------------------- */
-template <typename T, UInt ndim, class RetType> class TensorStorage {
+template <typename T, UInt ndim, class RetType>
+class TensorStorage : public TensorTrait {
public:
using value_type = T;
+ friend class Array<T>;
protected:
template <class TensorType> 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<T, ndim, RetType> & proxy) {
this->copySize(proxy);
this->values = proxy.storage();
this->wrapped = true;
}
protected:
TensorStorage(const TensorStorage & src) = delete;
public:
TensorStorage(const TensorStorage & src, bool deep_copy)
: values(nullptr), wrapped(false) {
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<T>{},
+ "Cannot create a tensor on non trivial types");
DimHelper<ndim>::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<ndim>::setDims(m, n, p, this->n);
this->computeSize();
this->values = data;
this->wrapped = true;
}
public:
/* ------------------------------------------------------------------------ */
template <class TensorType> inline void shallowCopy(const TensorType & src) {
this->copySize(src);
if (!this->wrapped)
delete[] this->values;
this->values = src.storage();
this->wrapped = true;
}
/* ------------------------------------------------------------------------ */
template <class TensorType> inline void deepCopy(const TensorType & src) {
this->copySize(src);
+
if (!this->wrapped)
delete[] this->values;
+
+ static_assert(std::is_trivially_constructible<T>{},
+ "Cannot create a tensor on non trivial types");
this->values = new T[this->_size];
+
+ static_assert(std::is_trivially_copyable<T>{},
+ "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<RetType &>(src));
}
/* ------------------------------------------------------------------------ */
inline TensorStorage & operator=(const RetType & src) {
if (this != &src) {
if (this->wrapped) {
+ static_assert(std::is_trivially_copyable<T>{},
+ "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 <class R>
inline RetType & operator+=(const TensorStorage<T, ndim, R> & 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<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
template <class R>
inline RetType & operator-=(const TensorStorage<T, ndim, R> & 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<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator+=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i)
*(a++) += x;
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator-=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i)
*(a++) -= x;
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator*=(const T & x) {
T * a = this->storage();
for (UInt i = 0; i < _size; ++i)
*(a++) *= x;
return *(static_cast<RetType *>(this));
}
/* ---------------------------------------------------------------------- */
inline RetType & operator/=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i)
*(a++) /= x;
return *(static_cast<RetType *>(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<RetType *>(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 <class TensorType> 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:
- friend class Array<T>;
-
inline void computeSize() {
_size = 1;
for (UInt d = 0; d < ndim; ++d)
_size *= this->n[d];
}
protected:
template <typename R, NormType norm_type> struct NormHelper {
template <class Ten> 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 <typename R> struct NormHelper<R, L_1> {
template <class Ten> 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 <typename R> struct NormHelper<R, L_2> {
template <class Ten> 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 <typename R> struct NormHelper<R, L_inf> {
template <class Ten> 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<L_p> @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 <NormType norm_type> inline T norm() const {
return NormHelper<T, norm_type>::norm(*this);
}
protected:
UInt n[ndim];
UInt _size;
T * values;
bool wrapped{false};
};
-// template <typename T, UInt ndim, class RetType>
-// inline TensorProxy<T, ndim, RetType>::TensorProxy(
-// const TensorStorage<T, ndim, RetType> & other) {
-// this->values = other.storage();
-// for (UInt i = 0; i < ndim; ++i)
-// this->n[i] = other.size(i);
-// }
-
/* -------------------------------------------------------------------------- */
/* Vector */
/* -------------------------------------------------------------------------- */
template <typename T> class Vector : public TensorStorage<T, 1, Vector<T>> {
typedef TensorStorage<T, 1, Vector<T>> parent;
public:
using value_type = typename parent::value_type;
using proxy = VectorProxy<T>;
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 VectorProxy<T> & src) : parent(src) {}
+ Vector(const TensorProxy<T, 1, Vector> & src) : parent(src) {}
Vector(std::initializer_list<T> list) : parent(list.size(), 0, 0, T()) {
UInt i = 0;
for (auto val : list) {
operator()(i++) = val;
}
}
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<T> & operator*=(Real x) { return parent::operator*=(x); }
inline Vector<T> & operator/=(Real x) { return parent::operator/=(x); }
/* ------------------------------------------------------------------------ */
inline Vector<T> & operator*=(const Vector<T> & vect) {
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<T> & 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<T> & v1, const Vector<T> & 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 void solve(const Matrix<T> & A, const Vector<T> & 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 <bool tr_A>
inline void mul(const Matrix<T> & A, const Vector<T> & x, Real alpha = 1.0);
/* ------------------------------------------------------------------------ */
inline Real norm() const { return parent::template norm<L_2>(); }
template <NormType nt> inline Real norm() const {
return parent::template norm<nt>();
}
/* ------------------------------------------------------------------------ */
inline void normalize() {
Real n = norm();
operator/=(n);
}
/* ------------------------------------------------------------------------ */
/// norm of (*this - x)
inline Real distance(const Vector<T> & 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<T> & 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<T> & 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<T> & v) const { return equal(v); }
inline bool operator!=(const Vector<T> & v) const { return !operator==(v); }
inline bool operator<(const Vector<T> & v) const { return compare(v) == -1; }
inline bool operator>(const Vector<T> & v) const { return compare(v) == 1; }
/* ------------------------------------------------------------------------ */
/// 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 << "]";
}
// friend class ::akantu::Array<T>;
};
using RVector = Vector<Real>;
/* ------------------------------------------------------------------------ */
template <>
inline bool Vector<UInt>::equal(const Vector<UInt> & 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;
}
/* ------------------------------------------------------------------------ */
/* Matrix */
/* ------------------------------------------------------------------------ */
template <typename T> class Matrix : public TensorStorage<T, 2, Matrix<T>> {
typedef TensorStorage<T, 2, Matrix<T>> parent;
public:
using value_type = typename parent::value_type;
using proxy = MatrixProxy<T>;
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<T> & src) : parent(src) {}
Matrix(std::initializer_list<std::initializer_list<T>> list) {
+ static_assert(std::is_trivially_copyable<T>{},
+ "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;
}
}
virtual ~Matrix() = 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<T> 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<T>(this->values + j * this->n[0], this->n[0]);
}
inline const VectorProxy<T> 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<T>(this->values + j * this->n[0], this->n[0]);
}
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<false, false>(*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<false, false>(C, B);
return *this;
}
/* ---------------------------------------------------------------------- */
template <bool tr_A, bool tr_B>
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<tr_A, tr_B>(this->rows(), this->cols(), k, alpha, A.storage(),
B.storage(), 0., this->storage());
}
/* ---------------------------------------------------------------------- */
inline void outerProduct(const Vector<T> & A, const Vector<T> & 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<T> & eigs) : eigs(eigs) {}
bool operator()(const UInt & a, const UInt & b) const {
return (eigs(a) > eigs(b));
}
private:
const Vector<T> & eigs;
};
public:
/* ---------------------------------------------------------------------- */
inline void eig(Vector<T> & eigenvalues, Matrix<T> & 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() != NULL)
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<T> tmp = *this;
Vector<T> tmp_eigs(eigenvalues.size());
Matrix<T> 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<UInt> 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<T> & eigenvalues) const {
Matrix<T> 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<T> eye(UInt m, T alpha = 1.) {
Matrix<T> 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 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<T> & 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_DEBUG_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 <typename T>
template <bool tr_A>
inline void Vector<T>::mul(const Matrix<T> & A, const Vector<T> & 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<tr_A>(A.rows(), A.cols(), alpha, A.storage(), x.storage(),
0., this->storage());
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline std::ostream & operator<<(std::ostream & stream,
const Matrix<T> & _this) {
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline std::ostream & operator<<(std::ostream & stream,
const Vector<T> & _this) {
_this.printself(stream);
return stream;
}
/* ------------------------------------------------------------------------ */
/* Tensor3 */
/* ------------------------------------------------------------------------ */
template <typename T> class Tensor3 : public TensorStorage<T, 3, Tensor3<T>> {
typedef TensorStorage<T, 3, Tensor3<T>> parent;
public:
using value_type = typename parent::value_type;
using proxy = Tensor3Proxy<T>;
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) {}
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<T> 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<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline const MatrixProxy<T> 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<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline MatrixProxy<T> operator[](UInt k) {
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline const MatrixProxy<T> operator[](UInt k) const {
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
};
/* -------------------------------------------------------------------------- */
// support operations for the creation of other vectors
/* -------------------------------------------------------------------------- */
template <typename T>
Vector<T> operator*(const T & scalar, const Vector<T> & a) {
Vector<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Vector<T> operator*(const Vector<T> & a, const T & scalar) {
Vector<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Vector<T> operator/(const Vector<T> & a, const T & scalar) {
Vector<T> r(a);
r /= scalar;
return r;
}
template <typename T>
Vector<T> operator*(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r *= b;
return r;
}
template <typename T>
Vector<T> operator+(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r += b;
return r;
}
template <typename T>
Vector<T> operator-(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r -= b;
return r;
}
template <typename T>
Vector<T> operator*(const Matrix<T> & A, const Vector<T> & b) {
Vector<T> r(b.size());
r.template mul<false>(A, b);
return r;
}
/* -------------------------------------------------------------------------- */
template <typename T>
Matrix<T> operator*(const T & scalar, const Matrix<T> & a) {
Matrix<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Matrix<T> operator*(const Matrix<T> & a, const T & scalar) {
Matrix<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Matrix<T> operator/(const Matrix<T> & a, const T & scalar) {
Matrix<T> r(a);
r /= scalar;
return r;
}
template <typename T>
Matrix<T> operator+(const Matrix<T> & a, const Matrix<T> & b) {
Matrix<T> r(a);
r += b;
return r;
}
template <typename T>
Matrix<T> operator-(const Matrix<T> & a, const Matrix<T> & b) {
Matrix<T> r(a);
r -= b;
return r;
}
-} // akantu
+} // namespace akantu
#endif /* __AKANTU_AKA_TYPES_HH__ */