Diffusion akantu (bugfixes/27-bug-in-the-python-interface) Browse

27-bug-in-the…

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: Tue Aug 19 2014
			*
			* @brief description of the "simple" types
			*
			* @section LICENSE
			*
			* Copyright (©) 2010-2012, 2014 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 "aka_array.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)
			*/
			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 */
			/* -------------------------------------------------------------------------- */

			/**
			* 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
			*/
			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);

			public:
			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 TensorProxy & 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 TensorProxy & operator=(const TensorProxy & 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;
			}

			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) { }
			VectorProxy & operator=(const type & src) {
			parent::operator=(src);
			return *this;
			}
			VectorProxy & operator=(const VectorProxy & src) {
			parent::operator=(src);
			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> > {
			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) { }
			MatrixProxy & operator=(const type & src) {
			parent::operator=(src);
			return *this;
			}
			MatrixProxy & operator=(const MatrixProxy & src) {
			parent::operator=(src);
			return *this;
			}
			};

			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) { }
			Tensor3Proxy & operator=(const type & src) {
			parent::operator=(src);
			return *this;
			}
			Tensor3Proxy & operator=(const Tensor3Proxy & src) {
			parent::operator=(src);
			return *this;
			}
			};

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

			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(this->values, 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(this->values, 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));
			}

			/* ------------------------------------------------------------------------ */
			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->values + i); }
			inline const T& operator[](UInt i) const { return *(this->values + 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 + jthis->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 + jthis->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 + (kthis->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 + (kthis->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 + kthis->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 + kthis->n[0]this->n[1], this->n[0], this->n[1]);
			}

			inline MatrixProxy<T> operator[](UInt k) {
			return Matrix<T>(this->values + kthis->n[0]this->n[1], this->n[0], this->n[1]);
			}

			inline const MatrixProxy<T> operator[](UInt k) const {
			return MatrixProxy<T>(this->values + kthis->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__ */

rAKA/src/common/aka_types.hhf58b7e7a1b1ebugfixes/27-bug-in-the…

aka_types.hhView Options

rAKA/src/common/aka_types.hh
f58b7e7a1b1ebugfixes/27-bug-in-the…

aka_types.hh
View Options