aka_types.hh
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, May 2, 05:49

aka_types.hh
View Options

	/**
	* @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 <iomanip>
	/* -------------------------------------------------------------------------- */

	#ifndef __AKANTU_AKA_TYPES_HH__
	#define __AKANTU_AKA_TYPES_HH__

	__BEGIN_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 */
	/* -------------------------------------------------------------------------- */

	/**
	* @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:
	TensorProxy(T * data, UInt m, UInt n, UInt p) {
	DimHelper<ndim>::setDims(m, n, p, this->n);
	this->values = data;
	}

	TensorProxy(const TensorProxy & other) {
	this->values = other.storage();
	for (UInt i = 0; i < ndim; ++i)
	this->n[i] = other.n[i];
	}

	inline TensorProxy(const TensorStorage<T, ndim, RetType> & other);

	typedef typename RetType::proxy RetTypeProxy;

	public:
	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; }

	inline RetTypeProxy & operator=(const RetType & src) {
	AKANTU_DEBUG_ASSERT(
	src.size() == this->size(),
	"You are trying to copy two tensors with different sizes");
	memcpy(this->values, src.storage(), this->size() * sizeof(T));
	return *this;
	}

	inline RetTypeProxy & operator=(const RetTypeProxy & src) {
	AKANTU_DEBUG_ASSERT(
	src.size() == this->size(),
	"You are trying to copy two tensors with different sizes");
	memcpy(this->values, src.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> > {
	typedef TensorProxy<T, 1, Vector<T> > parent;
	typedef Vector<T> type;

	public:
	VectorProxy(T * data, UInt n) : parent(data, n, 0, 0) {}
	VectorProxy(const VectorProxy & src) : parent(src) {}
	VectorProxy(const Vector<T> & src) : parent(src) {}
	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> > {
	typedef TensorProxy<T, 2, Matrix<T> > parent;
	typedef Matrix<T> type;

	public:
	MatrixProxy(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
	MatrixProxy(const MatrixProxy & src) : parent(src) {}
	MatrixProxy(const type & src) : parent(src) {}
	};

	template <typename T>
	class Tensor3Proxy : public TensorProxy<T, 3, Tensor3<T> > {
	typedef TensorProxy<T, 3, Tensor3<T> > parent;
	typedef Tensor3<T> type;

	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) {}
	};

	/* -------------------------------------------------------------------------- */
	/* Tensor base class */
	/* -------------------------------------------------------------------------- */
	template <typename T, UInt ndim, class RetType> class TensorStorage {
	public:
	typedef T value_type;

	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(NULL), wrapped(false) {
	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(NULL), wrapped(false) {
	if (deep_copy)
	this->deepCopy(src);
	else
	this->shallowCopy(src);
	}

	protected:
	TensorStorage(UInt m, UInt n, UInt p, const T & def) {
	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;
	this->values = new T[this->_size];
	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 RetType & src) {
	if (this != &src) {
	if (this->wrapped) {
	// 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;
	};

	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:
	typedef typename parent::value_type value_type;
	typedef VectorProxy<T> proxy;

	public:
	Vector() : parent() {}
	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) {}

	public:
	virtual ~Vector(){};

	/* ------------------------------------------------------------------------ */
	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(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 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>;
	};

	typedef Vector<Real> RVector;

	/* ------------------------------------------------------------------------ */
	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:
	typedef typename parent::value_type value_type;
	typedef MatrixProxy<T> proxy;

	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) {}

	virtual ~Matrix() {}
	/* ------------------------------------------------------------------------ */
	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:
	typedef typename parent::value_type value_type;
	typedef Tensor3Proxy<T> proxy;

	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>
	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;
	}

	__END_AKANTU__

	#endif /* __AKANTU_AKA_TYPES_HH__ */

aka_types.hhNo OneTemporaryActions

File Metadata

aka_types.hhView Options

Event Timeline

aka_types.hh
No OneTemporary
Actions

aka_types.hh
View Options