/****************************************************************************** * * AMDiS - Adaptive multidimensional simulations * * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved. * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis * * Authors: * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al. * * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. * * * This file is part of AMDiS * * See also license.opensource.txt in the distribution. * ******************************************************************************/ /** \file MatrixVectorOperations.h */ #ifndef AMDIS_MATVEC_OPERATIONS_H #define AMDIS_MATVEC_OPERATIONS_H #include "Traits.h" namespace AMDiS { // --------------------------------------------------------------------------- // Operations with Vector and Matrix /// Matrix vector multiplication. template inline const Vector& mv(const Matrix& m, const Vector& v, Vector& result) { TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n"); TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n"); T *resultIt; T const* mIt; T const* vIt; for (resultIt = result.begin(), mIt = m.begin(); resultIt != result.end(); ++resultIt) { *resultIt = 0; for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt) *resultIt += *mIt * *vIt; } return result; } /// Vector addition. template inline const Vector& add(const Vector& v1, const Vector& v2, Vector& result) { TEST_EXIT_DBG(v1.getSize() == v2.getSize())("invalid size in test v1 == v2\n"); TEST_EXIT_DBG(v2.getSize() == result.getSize())("invalid size in test v2 == result\n"); T const* v1It; T const* v2It; T* resultIt; for (v1It = v1.begin(), v2It = v2.begin(), resultIt = result.begin(); v1It != v1.end(); ++v1It, ++v2It, ++resultIt) *resultIt = *v1It + *v2It; return result; } /// scalar * vector template inline const Vector& mult(const S& scal, const Vector& v, Vector& result) { TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n"); T const* vIt; T* resultIt; for (vIt = v.begin(), resultIt = result.begin(); vIt != v.end(); ++vIt, ++resultIt) *resultIt = scal * *vIt; return result; } /// vector + scalar template inline const Vector& add(const Vector& v, const T& scal, Vector& result) { TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n"); T const* vIt; T* resultIt; for (vIt = v.begin(), resultIt = result.begin(); vIt != v.end(); ++vIt, ++resultIt) *resultIt = *vIt + scal; return result; } /// y = a * x + y. template inline const Vector& axpy(const S& a, const Vector &x, Vector &y) { TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n"); T const* xIt; T* yIt; for (xIt = x.begin(), yIt = y.begin(); xIt != x.end(); ++xIt, ++yIt) *yIt += a * *xIt; return y; } // times / divides // --------------- /// Matrix vector multiplication: vector := matrix * vector template inline const Vector& operator*=(const Vector& v, const Matrix& m) { return mv(m, v, v); } /// Matrix vector multiplication: vector := matrix * vector template inline Vector operator*(const Matrix& m, const Vector& v) { Vector result(m.getNumCols()); return mv(m, v, result); } /// Scalar product: scalar := vector * vector template inline typename traits::mult_type::type operator*(const Vector& v1, const Vector& v2) { typename traits::mult_type::type result; nullify(result); T const* v1It; S const* v2It; for (v1It = v1.begin(), v2It = v2.begin(); v1It != v1.end(); ++v1It, ++v2It) result += *v1It * *v2It; return result; } /// vector *= scalar (elementwise) template typename enable_if< traits::is_multiplicable, Vector >::type & operator*=(Vector& v, S scal) { T *vIt; for (vIt = v.begin(); vIt != v.end(); ++vIt) *vIt *= scal; return v; } /// vector := vector * scalar (elementwise) template typename enable_if< traits::is_multiplicable, Vector >::type operator*(Vector result, S scal) { result *= scal; return result; } /// vector /= scalar (elementwise) template typename enable_if< traits::is_multiplicable, Vector >::type & operator/=(Vector& v, S scal) { T *vIt; for (vIt = v.begin(); vIt != v.end(); ++vIt) *vIt /= scal; return v; } /// vector := vector / scalar (elementwise) template typename enable_if< traits::is_multiplicable, Vector >::type operator/(Vector result, S scal) { result /= scal; return result; } // plus / minus // ------------ /// vector += scalar template typename boost::enable_if::type, Vector >::type& operator+=(Vector& x, S scal) { T *xIt; for (xIt = x.begin(); xIt != x.end(); ++xIt) *xIt += scal; return x; } /// vector := vector + scalar template typename boost::enable_if::type, Vector >::type operator+(Vector result, T scal) { result += scal; return result; } /// vector += vector template Vector& operator+=(Vector& x, const Vector& y) { T *xIt; S const* yIt; for (xIt = x.begin(), yIt = y.begin(); xIt != x.end(); ++xIt, ++yIt) *xIt += *yIt; return x; } /// vector := vector + vector template Vector operator+(Vector result, const Vector& v2) { result += v2; return result; } /// vector -= vector template Vector& operator-=(Vector& x, const Vector& y) { T* xIt; S const* yIt; for (xIt = x.begin(), yIt = y.begin(); xIt != x.end(); ++xIt, ++yIt) *xIt -= *yIt; return x; } /// vector := vector - vector template Vector operator-(Vector result, const Vector& v2) { result -= v2; return result; } // special operators // ----------------- /// 2-norm of a vector template inline T norm(const Vector *v) { T result; nullify(result); for (T const* vIt = v->begin(); vIt != v->end(); ++vIt) result += *vIt * *vIt; return std::sqrt(result); } /// 2-norm of a vector template inline T norm(const Vector& v) { return norm(&v); } /// cross-product of two vectors: vector := vector x vector (only in 3d) template void vectorProduct(const Vector& x, const Vector& y, Vector& z) { FUNCNAME_DBG("vectorProduct()"); TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n"); z[0] = x[1] * y[2] - x[2] * y[1]; z[1] = x[2] * y[0] - x[0] * y[2]; z[2] = x[0] * y[1] - x[1] * y[0]; } // --------------------------------------------------------------------------- // Operations with WorldVector and WorldMatrix // times / divides // --------------- /// Scalar product: scalar := vector * vector template inline typename traits::mult_type::type operator*(const WorldVector& v1, const WorldVector& v2) { typename traits::mult_type::type result; nullify(result); T const* v1It; S const* v2It; for (v1It = v1.begin(), v2It = v2.begin(); v1It != v1.end(); ++v1It, ++v2It) result += *v1It * *v2It; return result; } /// vector := vector * scalar (elementwise) template typename enable_if< traits::is_multiplicable, WorldVector >::type operator*(WorldVector const& v, S scal) { WorldVector result = v; result *= scal; // calls operator*=(Vector, S) return result; } /// vector := scalar * vector (elementwise) template typename enable_if< traits::is_multiplicable, WorldVector >::type operator*(S scal, WorldVector const& v) { WorldVector result = v; result *= scal; // calls operator*=(Vector, S) return result; } /// vector := vector / scalar (elementwise) template typename enable_if< traits::is_multiplicable, WorldVector >::type operator/(WorldVector const& v, S scal) { WorldVector result = v; result /= scal; // calls operator/=(Vector, S) return result; } /// matrix *= scalar (elementwise) // template // WorldMatrix& operator*=(WorldMatrix& m, T scal) // { // for (T* mIt = m.begin(); mIt != m.end(); mIt++) // *mIt *= scal; // // return m; // } /// matrix := matrix * scalar (elementwise) template typename enable_if< traits::is_multiplicable, WorldMatrix >::type operator*(WorldMatrix const& m, S scal) { WorldMatrix result = m; result *= scal; // calls operator*=(Vector, S) return result; } /// matrix := scalar * matrix (elementwise) template typename enable_if< traits::is_multiplicable, WorldMatrix >::type operator*(S scal, WorldMatrix const& m) { WorldMatrix result = m; result *= scal; // calls operator*=(Vector, S) return result; } /// matrix := matrix / scalar (elementwise) template typename enable_if< traits::is_multiplicable, WorldMatrix >::type operator/(WorldMatrix const& m, S scal) { WorldMatrix result = m; result /= scal; // calls operator/=(Vector, S) return result; } /// vector := matrix * vector template WorldVector operator*(const WorldMatrix& M, const WorldVector& v ) { WorldVector res; res.multMatrixVec(M,v); return res; } /// matrix := matrix * matrix template WorldVector > operator*(const WorldVector >& A, const WorldVector >& B) { WorldVector > result; nullify(result); for (size_t r = 0; r < num_rows(A); r++) for (size_t c = 0; c < num_cols(A); c++) for (size_t i = 0; i < num_cols(A); i++) result[r][c] += A[r][i] * B[i][c]; return result; } // plus / minus // ------------ // NOTE: call operators of Vector directly #if 0 template WorldVector& operator+=(WorldVector& v1, const WorldVector& v2) { add(v1, v2, v1); return v1; } template WorldVector& operator-=(WorldVector& v1, const WorldVector& v2) { axpy(-1.0, v2, v1); return v1; } #endif /// vector := vector + vector template WorldVector& operator+=(WorldVector& v1, WorldVector const& v2) { static_cast&>(v1) += static_cast const&>(v2); return v1; } /// vector := vector + vector template WorldVector::type> operator+(WorldVector result, const WorldVector& v2) { result += v2; // calls operator+=(Vector, Vector) return result; } /// vector := vector - vector template WorldVector::type> operator-(WorldVector result, const WorldVector& v2) { result -= v2; // calls operator-=(Vector, Vector) return result; } /// matrix += matrix template WorldMatrix& operator+=(WorldMatrix& m1, const WorldMatrix& m2) { T* m1It; S const* m2It; for (m1It = m1.begin(), m2It = m2.begin(); m1It != m1.end(); m1It++, m2It++) *m1It += *m2It; return m1; } /// matrix += matrix template Matrix& operator+=(Matrix& m1, const Matrix& m2) { T* m1It; S const* m2It; for (m1It = m1.begin(), m2It = m2.begin(); m1It != m1.end(); m1It++, m2It++) *m1It += *m2It; return m1; } /// matrix := matrix + matrix template WorldMatrix operator+(WorldMatrix M1, const WorldMatrix& M2 ) { M1 += M2; return M1; } /// matrix := matrix + matrix template Matrix operator+(Matrix M1, const Matrix& M2 ) { M1 += M2; return M1; } /// matrix -= matrix template WorldMatrix& operator-=(WorldMatrix& m1, const WorldMatrix& m2) { T *m1It; S const* m2It; for (m1It = m1.begin(), m2It = m2.begin(); m1It != m1.end(); m1It++, m2It++) *m1It -= *m2It; return m1; } /// matrix -= matrix template Matrix& operator-=(Matrix& m1, const Matrix& m2) { T *m1It; S const* m2It; for (m1It = m1.begin(), m2It = m2.begin(); m1It != m1.end(); m1It++, m2It++) *m1It -= *m2It; return m1; } /// matrix := matrix - matrix template WorldMatrix operator-(WorldMatrix M1, const WorldMatrix& M2 ) { M1 -= M2; return M1; } /// matrix := matrix - matrix template Matrix operator-(Matrix M1, const Matrix& M2 ) { M1 -= M2; return M1; } // unary minus operators // --------------------- /// vector := -vector (elementwise) template WorldVector operator-(WorldVector v) { v *= -1.0; return v; } /// matrix := -matrix (elementwise) template WorldMatrix operator-(WorldMatrix v) { v *= -1.0; return v; } // comparison operators // -------------------- /// test for less-then (elementwise) up to DBL_TOL inline bool operator<(const WorldVector& v1, const WorldVector& v2) { int dow = Global::getGeo(WORLD); for (int i = 0; i < dow; i++) { if (std::abs(v1[i] - v2[i]) < DBL_TOL) continue; return v1[i] < v2[i]; } return false; } /// test for equality (elementwise) up to DBL_TOL inline bool operator==(const WorldVector& v1, const WorldVector& v2) { int dow = Global::getGeo(WORLD); for (int i = 0; i < dow; i++) if (std::abs(v1[i] - v2[i]) > DBL_TOL) return false; return true; } /// Comparison operator for Vector. template inline bool operator==(Vector const& lhs, Vector const& rhs) { if (lhs.getSize() != rhs.getSize()) return false; T const* rhsIt; T const* lhsIt; for (rhsIt = rhs.begin(), lhsIt = lhs.begin(); rhsIt != rhs.end(); ++rhsIt, ++lhsIt) if (*lhsIt != *rhsIt) return false; return true; } /// Comparison operator for Vector. template inline bool operator!=(Vector const& lhs, Vector const& rhs) { return !(lhs == rhs); } /// Comparison operator for Matrix. template inline bool operator==(Matrix const& lhs, Matrix const& rhs) { if (lhs.getNumRows() != rhs.getNumRows()) return false; if (lhs.getNumCols() != rhs.getNumCols()) return false; return (static_cast const&>(lhs) == static_cast const&>(rhs)); } /// Comparison operator for Matrix. template inline bool operator!=(Matrix const& lhs, Matrix const& rhs) { return !(lhs == rhs); } // special operators // ----------------- /// wrapper for nullify template void set_to_zero(WorldVector >& mat) { nullify(mat); } // NOTE: call norm(Vector) directly #if 0 inline double norm(const WorldVector& v) { double val = 0.0; for (int i = 0; i < Global::getGeo(WORLD); i++) val += v[i] * v[i]; return sqrt(val); } #endif /// returns the euclidian distance of a and b template double absteukl(const FixVec& a,const FixVec& b) { double erg = 0.0; for (int i = 0; i < a.getSize(); ++i) erg += sqr(a[i] - b[i]); return std::sqrt(erg); } } // end namespace AMDiS #endif // AMDIS_MATVEC_OPERATIONS_H