/** \file VectorOperations.h */ #ifndef VECTOR_OPERATIONS_H #define VECTOR_OPERATIONS_H #include "AMDiS.h" #include "ValueTypes.h" #if HAVE_PARALLEL_DOMAIN_AMDIS #include "parallel/StdMpi.h" #endif #include #include #include #include #include #include #include #include using namespace AMDiS; template const WorldVector& operator*=(WorldVector& v1, const WorldVector& v2) { T* v1It, *v2It; for (v1It = v1.begin(), v2It = v2.begin(); v1It != v1.end(); v1It++, v2It++) *v1It *= *v2It; return v1; } // x ~= y template class compareTol : public binary_function { public: compareTol(double tol_) : tol(tol_) {} bool operator()(const T &x, const T &y) const { return sqrt((x-y)*(x-y)) struct comparePair : public binary_function, std::pair, bool> { bool operator()(const std::pair &x, const std::pair &y) const { return (component==0?(x.first struct compare0 : public binary_function, std::pair, bool> { bool operator()(const std::pair &x, const std::pair &y) const { return x.first struct compare1 : public binary_function, std::pair, bool> { bool operator()(const std::pair &x, const std::pair &y) const { return x.second, double>, std::pair, double>, bool> { public: SortPair(int component_=1,int idx_=0) : component(component_),idx(idx_) {} bool operator()(const std::pair, double> &x, const std::pair, double> &y) const { return (component==1?(x.first)[idx]<(y.first)[idx]:x.second struct num_rows {}; /// size implementation for STL vectors template struct num_rows< std::vector > { typedef std::size_t type; type operator()(const std::vector& v) { return v.size(); } }; /// size implementation for (1D) arrays interpreted as vectors template struct num_rows { typedef std::size_t type; type operator()(const Value[Size]) { return Size; } }; /// size implementation for (2D and higher) arrays interpreted as matrices template struct num_rows { typedef std::size_t type; type operator()(const Value[Rows][Cols]) { return Rows; } }; /// size implementation for AMDiS WorldVectors template struct num_rows< WorldVector > { typedef std::size_t type; type operator()(const WorldVector& v) { return static_cast(v.getSize()); } }; /// size implementation for AMDiS WorldMatrices template struct num_rows< WorldMatrix > { typedef std::size_t type; type operator()(const WorldMatrix& v) { return static_cast(v.getNumRows()); } }; /// size implementation for arithmetic types (double, float, int, ...) template struct num_rows< Value, typename boost::enable_if< boost::is_arithmetic >::type > { typedef std::size_t type; type operator()(const Value& v) { return 1; } }; //________________________________________________________________________________ /// General declaration, used to disable unsupported types template struct num_cols {}; /// size implementation for AMDiS WorldMatrices template struct num_cols< WorldMatrix > { typedef std::size_t type; type operator()(const WorldMatrix& v) { return static_cast(v.getNumCols()); } }; //________________________________________________________________________________ /// General declaration, used to disable unsupported types template struct resize {}; /// change_dim implementation for AMDiS WorldVectors template struct resize< WorldVector > { typedef void vector_void_type; void operator()(WorldVector& v, size_t r) { TEST_EXIT(Global::getGeo(WORLD) == r) ("WorldVectors can not be resized!\n"); } }; /// change_dim implementation for STL vectors template struct resize< std::vector > { typedef void vector_void_type; void operator()(std::vector& v, size_t r) { v.resize(r); } }; /// change_dim implementation for MTL4 vectors template struct resize< mtl::dense_vector > { typedef void vector_void_type; void operator()(mtl::dense_vector& v, size_t r) { v.change_dim(r); } }; // _________________________________________________________________________________ /// change_dim implementation for AMDiS WorldMatrices template struct resize< WorldMatrix > { typedef void matrix_void_type; void operator()(WorldMatrix& v, size_t r, size_t c) { size_t dow = static_cast(Global::getGeo(WORLD)); TEST_EXIT(dow == r && dow == c) ("WorldMatrices can not be resized!\n"); } }; /// change_dim implementation for MTL4 matrix template struct resize< mtl::matrix::base_matrix > { typedef void matrix_void_type; void operator()(mtl::matrix::base_matrix& v, size_t r, size_t c) { v.change_dim(r,c); } }; /// change_dim implementation for MTL4 matrix template struct resize< mtl::matrix::dense2D > { typedef void matrix_void_type; void operator()(mtl::matrix::dense2D& v, size_t r, size_t c) { v.change_dim(r,c); } }; /// change_dim implementation for MTL4 matrix template struct resize< mtl::matrix::compressed2D > { typedef void matrix_void_type; void operator()(mtl::matrix::compressed2D& v, size_t r, size_t c) { v.change_dim(r,c); } }; // _________________________________________________________________________________ /// General declaration, used to disable unsupported types template struct at {}; template struct at< WorldMatrix > { typedef Value type; type& operator()(WorldMatrix& v, size_t r, size_t c) { return v[r][c]; } type const& operator()(const WorldMatrix& v, size_t r, size_t c) { return v[r][c]; } }; template struct at< mtl::matrix::base_matrix > { typedef Value type; type& operator()(mtl::matrix::base_matrix& v, size_t r, size_t c) { return v(r,c); } type const& operator()(const mtl::matrix::base_matrix& v, size_t r, size_t c) { return v(r,c); } }; template struct at< mtl::matrix::dense2D > { typedef Value type; type& operator()(mtl::matrix::dense2D& v, size_t r, size_t c) { return v(r,c); } type const& operator()(const mtl::matrix::dense2D& v, size_t r, size_t c) { return v(r,c); } }; template struct at< mtl::matrix::compressed2D > { typedef Value type; type& operator()(mtl::matrix::compressed2D& v, size_t r, size_t c) { return v(r,c); } type const& operator()(const mtl::matrix::compressed2D& v, size_t r, size_t c) { return v(r,c); } }; } /// num_rows function for non-MTL types (uses implicit enable_if) template typename traits::num_rows::type num_rows(const Collection& c) { return traits::num_rows()(c); } /// num_cols function for non-MTL types (uses implicit enable_if) template typename traits::num_cols::type num_cols(const Collection& c) { return traits::num_cols()(c); } /// resize function for vectors template typename traits::resize::vector_void_type resize(Collection& c, size_t rows) { traits::resize()(c, rows); } /// resize function for matrices template typename traits::resize::matrix_void_type resize(Collection& c, size_t rows, size_t cols) { traits::resize()(c, rows, cols); } /// at function for matrices template typename traits::at::type const& at(const Collection& c, size_t rows, size_t cols) { return traits::at()(c, rows, cols); } /// at function for matrices template typename traits::at::type& at(Collection& c, size_t rows, size_t cols) { return traits::at()(c, rows, cols); } /// 2-norm for vectors template typename boost::enable_if< is_vector, typename ValueType::type >::type norm(const Vector &b) { typename ValueType::type value; nullify(value); for (size_t i = 0; i < num_rows(b); ++i) value += sqr(b[i]); return sqrt(value); } /// 2-norm for MTL4 vectors template T norm(const mtl::dense_vector &b) { return two_norm(b); } /// 1-norm for matrices template typename boost::enable_if< is_matrix, double >::type matrix_norm(const Matrix &M) { typename Matrix::value_type asum; asum = 0.0; for (size_t i = 0; i < num_rows(M); ++i) { for (size_t j = 0; j < num_cols(M); ++j) asum += abs(M[i][j]); } return static_cast(asum); } /// 2-norm for vectors with eps-pertubation template typename boost::enable_if< is_vector, double >::type norm(const Vector &b, const double eps) { return std::max(static_cast(norm(b)), eps); } /// trace of a matrix template typename boost::enable_if< is_matrix, double >::type trace(Matrix &M) { TEST_EXIT_DBG(num_rows(M) == num_cols(M)) ("Matrix dimensions must be equal!\n"); typename Matrix::value_type sum; nullify(sum); for (size_t i = 0; i < num_rows(M); ++i) sum += M[i][i]; return static_cast(sum); } /// equivalent to matlab .* for vectors template typename boost::enable_if< boost::mpl::and_< is_vector, is_vector >, ResultVector >::type mult(const Vector1& v1, const Vector2& v2) { TEST_EXIT_DBG(num_rows(v1)==num_rows(v2)) ("Vectors must have the same size!\n"); ResultVector result(num_rows(v1)); for (size_t i = 0; i < num_rows(v1); ++i) result[i] = v1[i] * v2[i]; return result; } template typename boost::enable_if< boost::mpl::and_< is_matrix, is_vector >, ResultVector >::type mv(const Matrix &A, const Vector &x) { TEST_EXIT_DBG(num_cols(A)==num_rows(x)) ("Matrix and vector dimension must fit together!"); ResultVector result(num_rows(A)); for (size_t i = 0; i < num_rows(A); ++i) { result[i] = 0.0; for (size_t j = 0; j < num_cols(A); ++j) result[i] += at(A,i,j) * x[j]; } return result; } template typename boost::enable_if< boost::mpl::and_< is_matrix, is_vector, is_vector >, ResultVector >::type Axpy(const Matrix &A, const Vector1 &x, const Vector2 &y) { TEST_EXIT_DBG(num_cols(A)==num_row(x) && num_rows(A)==num_rows(y)) ("Matrix and vector dimensions must fit together!"); ResultVector result; resize(result, num_rows(A)); for (size_t i = 0; i < num_rows(A); ++i) { result[i] = y[i]; for (size_t j = 0; j < num_cols(A); ++j) result[i] += at(A,i,j) * x[j]; } return result; } // copy v -> w template typename boost::enable_if< boost::mpl::and_< is_vector, is_vector >, void >::type copy(const VectorIn &v, VectorOut &w) { for (size_t i = 0; i < std::min(num_rows(v), num_rows(w)); i++) w[i] = v[i]; } template typename boost::enable_if< boost::mpl::and_< is_matrix, is_matrix >, void >::type fillMatrix(const MatrixIn &m_in, MatrixOut &m) { resize(m, num_rows(m_in), num_cols(m_in)); for (size_t i = 0; i < num_rows(m_in); ++i) { for (size_t j = 0; j < num_cols(m_in); ++j) { at(m,i,j) = at(m_in,i,j); } } } template typename boost::enable_if< boost::mpl::and_< is_matrix, is_matrix >, void >::type copy(const MatrixIn &m_in, MatrixOut &m) { fillMatrix(m_in, m); } // v3:=[v1, v2] template typename boost::enable_if< boost::mpl::and_< is_vector, is_vector, is_vector >, void >::type merge(Vector1 &v1, Vector2 &v2, VectorOut &v3) { resize(v3, num_rows(v1) + num_rows(v2)); for (size_t i = 0; i < num_rows(v1); i++) v3[i] = v1[i]; for (size_t j = 0; j < num_rows(v2); j++) v3[j + num_rows(v1)] = v2[j]; } template void getMin(const Vector &v, typename ValueType::type &minVal, size_t &minIdx) { typedef typename ValueType::type T; TEST_EXIT(num_rows(v) > 0)("getMin of empty vector!\n"); minVal = v[0]; minIdx = 0; for (size_t i = 1; i < num_rows(v); i++) { if (v[i] < minVal) { minVal = v[i]; minIdx = i; } } } template void getMax(const Vector &v, typename ValueType::type &maxVal, size_t &maxIdx) { typedef typename ValueType::type T; TEST_EXIT(num_rows(v) > 0)("getMax of empty vector!\n"); maxVal = v[0]; maxIdx = 0; for (size_t i = 1; i < num_rows(v); i++) { if (v[i] > maxVal) { maxVal = v[i]; maxIdx = i; } } } template typename ValueType::type min(const Vector &vec) { typename ValueType::type minVal; size_t minIdx; getMin(vec, minVal, minIdx); return minVal; } template typename ValueType::type max(const Vector &vec) { typename ValueType::type maxVal; size_t maxIdx; getMax(vec, maxVal, maxIdx); return maxVal; } template typename ValueType::type absmin(const Vector &v) { typedef typename ValueType::type T; TEST_EXIT(num_rows(v) > 0)("absmin of empty vector!\n"); T minVal = abs(v[0]); for(size_t i = 1; i < num_rows(v); i++) { if (std::abs(v[i]) < minVal) { minVal = std::abs(v[i]); } } return minVal; } template typename ValueType::type absmax(const Vector &v) { typedef typename ValueType::type T; TEST_EXIT(num_rows(v) > 0)("absmax of empty vector!\n"); T maxVal = abs(v[0]); for(size_t i = 1; i < num_rows(v); i++) { if (std::abs(v[i]) > maxVal) { maxVal = std::abs(v[i]); } } return maxVal; } template typename ValueType::type sum(const Vector &v) { typedef typename ValueType::type T; size_t N = num_rows(v); T value; nullify(value); if (N > 0) { value = v[0]; for (size_t i = 1; i < N; ++i) value += v[i]; } return value; } template inline typename ValueType::type mean(const Vector &v) { double N = static_cast(num_rows(v)); return sum(v) * (1.0 / N); } template void getMean(const Vector &v, typename ValueType::type &meanValue) { meanValue = mean(v); } template void sort(std::vector &vec1, std::vector &vec2, Compare &comp) { std::vector > order(vec1.size()); // create vector of pairs for (size_t i=0; i void sort(std::vector > &vec, Compare &comp) { std::sort(vec.begin(), vec.end(), comp); } inline void unique(std::vector > &vec, double tol, std::vector &ind) { compareTol > comp(tol); unsigned newVec=0; for(unsigned i=0; i void unique(std::vector > &vec, double tol, int uniqueBy=1) { // comparePair comp = comparePair(tol,uniqueBy); // unsigned newVec=0; // for(unsigned i=0; i inline void filter(std::vector &x, const std::vector &filter) { TEST_EXIT_DBG(num_rows(x) == num_rows(filter)) ("number of rows of x and filter must be equal!\n"); size_t j = 0; for (size_t i = 0; i < num_rows(x); ++i) { if (filter[i]) { x[j] = x[i]; j++; } } x.erase(x.begin()+j,x.end()); } template void filter(std::vector &x, const std::vector &filter) { std::vector x_(num_rows(filter)); for (size_t i = 0; i < num_rows(filter); ++i) x_[i] = x[filter[i]]; swap(x, x_); } template typename ProductType::type, typename ValueType::type>::type dot(const Vector1 &vec1, const Vector2 &vec2) { typedef typename ProductType::type, typename ValueType::type>::type value_type; value_type value; nullify(value); size_t N = num_rows(vec1); if (N > 0) { value = vec1[0] * vec2[0]; for (size_t i = 1; i < N; ++i) value += vec1[i] * vec2[i]; } return value; } } #endif // VECTOR_OPERATIONS_H