From bc6d78982fb2b8d07827246e682acdf47d0e8944 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 18 Nov 2011 11:10:27 -0500 Subject: [PATCH] Bugs 157 and 377 - General tightening/testing of vectorwise ops: * add lots of static assertions making it very explicit when all these ops are supposed to work: ** all ops require the rhs vector to go in the right direction ** all ops already require that the lhs and rhs are of the same kind (matrix vs vector) otherwise we'd have to do complex work ** multiplicative ops (introduced Kibeom's patch) are restricted to arrays, if only because for matrices they could be ambiguous. * add a new test, vectorwiseop.cpp. * these compound-assign operators used to be implemented with for loops: for(Index j=0; j class VectorwiseOp typename ExtendedType::Type extendedTo(const DenseBase& other) const { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); + EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(Direction==Vertical, OtherDerived::MaxColsAtCompileTime==1), + YOU_PASSED_A_ROW_VECTOR_BUT_A_COLUMN_VECTOR_WAS_EXPECTED) + EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(Direction==Horizontal, OtherDerived::MaxRowsAtCompileTime==1), + YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED) return typename ExtendedType::Type (other.derived(), Direction==Vertical ? 1 : m_matrix.rows(), @@ -418,10 +421,9 @@ template class VectorwiseOp ExpressionType& operator=(const DenseBase& other) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) //eigen_assert((m_matrix.isNull()) == (other.isNull())); FIXME - for(Index j=0; j(m_matrix); + return const_cast(m_matrix = extendedTo(other.derived())); } /** Adds the vector \a other to each subvector of \c *this */ @@ -429,9 +431,8 @@ template class VectorwiseOp ExpressionType& operator+=(const DenseBase& other) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) - for(Index j=0; j(m_matrix); + EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) + return const_cast(m_matrix += extendedTo(other.derived())); } /** Substracts the vector \a other to each subvector of \c *this */ @@ -439,28 +440,29 @@ template class VectorwiseOp ExpressionType& operator-=(const DenseBase& other) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) - for(Index j=0; j(m_matrix); + EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) + return const_cast(m_matrix -= extendedTo(other.derived())); } - /** Multiplies the vector \a other to each subvector of \c *this */ + /** Multiples each subvector of \c *this by the vector \a other */ template ExpressionType& operator*=(const DenseBase& other) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) - for(Index j=0; j(m_matrix); } - /** Divides the vector \a other to each subvector of \c *this */ + /** Divides each subvector of \c *this by the vector \a other */ template ExpressionType& operator/=(const DenseBase& other) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) - for(Index j=0; j(m_matrix); } @@ -471,7 +473,8 @@ template class VectorwiseOp const typename ExtendedType::Type> operator+(const DenseBase& other) const { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) return m_matrix + extendedTo(other.derived()); } @@ -482,29 +485,36 @@ template class VectorwiseOp const typename ExtendedType::Type> operator-(const DenseBase& other) const { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) return m_matrix - extendedTo(other.derived()); } - /** Returns the expression of the multiplication of the vector \a other to each subvector of \c *this */ + /** Returns the expression where each subvector is the product of the vector \a other + * by the corresponding subvector of \c *this */ template EIGEN_STRONG_INLINE CwiseBinaryOp, const ExpressionTypeNestedCleaned, const typename ExtendedType::Type> operator*(const DenseBase& other) const { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType) + EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) return m_matrix * extendedTo(other.derived()); } - /** Returns the expression of the division between each subvector of \c *this and the vector \a other */ + /** Returns the expression where each subvector is the quotient of the corresponding + * subvector of \c *this by the vector \a other */ template CwiseBinaryOp, const ExpressionTypeNestedCleaned, const typename ExtendedType::Type> operator/(const DenseBase& other) const { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType) + EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) return m_matrix / extendedTo(other.derived()); } diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index d14c87c7c..672c6ae96 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h @@ -97,7 +97,10 @@ THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY, YOU_ARE_TRYING_TO_USE_AN_INDEX_BASED_ACCESSOR_ON_AN_EXPRESSION_THAT_DOES_NOT_SUPPORT_THAT, THIS_METHOD_IS_ONLY_FOR_1x1_EXPRESSIONS, - THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_OF_BOOL + THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_OF_BOOL, + THIS_METHOD_IS_ONLY_FOR_ARRAYS_NOT_MATRICES, + YOU_PASSED_A_ROW_VECTOR_BUT_A_COLUMN_VECTOR_WAS_EXPECTED, + YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED }; }; @@ -197,4 +200,15 @@ EIGEN_STATIC_ASSERT(internal::is_lvalue::value, \ THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY) +#define EIGEN_STATIC_ASSERT_ARRAYXPR(Derived) \ + EIGEN_STATIC_ASSERT((internal::is_same::XprKind, ArrayXpr>::value), \ + THIS_METHOD_IS_ONLY_FOR_ARRAYS_NOT_MATRICES) + +#define EIGEN_STATIC_ASSERT_SAME_XPR_KIND(Derived1, Derived2) \ + EIGEN_STATIC_ASSERT((internal::is_same::XprKind, \ + typename internal::traits::XprKind \ + >::value), \ + YOU_CANNOT_MIX_ARRAYS_AND_MATRICES) + + #endif // EIGEN_STATIC_ASSERT_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7086f5251..e1379d78b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -178,7 +178,7 @@ ei_add_test(dontalign) ei_add_test(evaluators) ei_add_test(sizeoverflow) ei_add_test(prec_inverse_4x4) - +ei_add_test(vectorwiseop) ei_add_test(simplicial_cholesky) ei_add_test(conjugate_gradient) diff --git a/test/vectorwiseop.cpp b/test/vectorwiseop.cpp new file mode 100644 index 000000000..d3518b7ec --- /dev/null +++ b/test/vectorwiseop.cpp @@ -0,0 +1,187 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Benoit Jacob +// +// Eigen 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. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen 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 or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#define EIGEN_NO_STATIC_ASSERT + +#include "main.h" + +template void vectorwiseop_array(const ArrayType& m) +{ + typedef typename ArrayType::Index Index; + typedef typename ArrayType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Array ColVectorType; + typedef Array RowVectorType; + + Index rows = m.rows(); + Index cols = m.cols(); + Index r = internal::random(0, rows-1), + c = internal::random(0, cols-1); + + ArrayType m1 = ArrayType::Random(rows, cols), + m2(rows, cols), + m3(rows, cols); + + ColVectorType colvec = ColVectorType::Random(rows); + RowVectorType rowvec = RowVectorType::Random(cols); + + // test addition + + m2 = m1; + m2.colwise() += colvec; + VERIFY_IS_APPROX(m2, m1.colwise() + colvec); + VERIFY_IS_APPROX(m2.col(c), m1.col(c) + colvec); + + VERIFY_RAISES_ASSERT(m2.colwise() += colvec.transpose()); + VERIFY_RAISES_ASSERT(m1.colwise() + colvec.transpose()); + + m2 = m1; + m2.rowwise() += rowvec; + VERIFY_IS_APPROX(m2, m1.rowwise() + rowvec); + VERIFY_IS_APPROX(m2.row(r), m1.row(r) + rowvec); + + VERIFY_RAISES_ASSERT(m2.rowwise() += rowvec.transpose()); + VERIFY_RAISES_ASSERT(m1.rowwise() + rowvec.transpose()); + + // test substraction + + m2 = m1; + m2.colwise() -= colvec; + VERIFY_IS_APPROX(m2, m1.colwise() - colvec); + VERIFY_IS_APPROX(m2.col(c), m1.col(c) - colvec); + + VERIFY_RAISES_ASSERT(m2.colwise() -= colvec.transpose()); + VERIFY_RAISES_ASSERT(m1.colwise() - colvec.transpose()); + + m2 = m1; + m2.rowwise() -= rowvec; + VERIFY_IS_APPROX(m2, m1.rowwise() - rowvec); + VERIFY_IS_APPROX(m2.row(r), m1.row(r) - rowvec); + + VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose()); + VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose()); + + // test multiplication + + m2 = m1; + m2.colwise() *= colvec; + VERIFY_IS_APPROX(m2, m1.colwise() * colvec); + VERIFY_IS_APPROX(m2.col(c), m1.col(c) * colvec); + + VERIFY_RAISES_ASSERT(m2.colwise() *= colvec.transpose()); + VERIFY_RAISES_ASSERT(m1.colwise() * colvec.transpose()); + + m2 = m1; + m2.rowwise() *= rowvec; + VERIFY_IS_APPROX(m2, m1.rowwise() * rowvec); + VERIFY_IS_APPROX(m2.row(r), m1.row(r) * rowvec); + + VERIFY_RAISES_ASSERT(m2.rowwise() *= rowvec.transpose()); + VERIFY_RAISES_ASSERT(m1.rowwise() * rowvec.transpose()); + + // test quotient + + m2 = m1; + m2.colwise() /= colvec; + VERIFY_IS_APPROX(m2, m1.colwise() / colvec); + VERIFY_IS_APPROX(m2.col(c), m1.col(c) / colvec); + + VERIFY_RAISES_ASSERT(m2.colwise() /= colvec.transpose()); + VERIFY_RAISES_ASSERT(m1.colwise() / colvec.transpose()); + + m2 = m1; + m2.rowwise() /= rowvec; + VERIFY_IS_APPROX(m2, m1.rowwise() / rowvec); + VERIFY_IS_APPROX(m2.row(r), m1.row(r) / rowvec); + + VERIFY_RAISES_ASSERT(m2.rowwise() /= rowvec.transpose()); + VERIFY_RAISES_ASSERT(m1.rowwise() / rowvec.transpose()); +} + +template void vectorwiseop_matrix(const MatrixType& m) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix ColVectorType; + typedef Matrix RowVectorType; + + Index rows = m.rows(); + Index cols = m.cols(); + Index r = internal::random(0, rows-1), + c = internal::random(0, cols-1); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2(rows, cols), + m3(rows, cols); + + ColVectorType colvec = ColVectorType::Random(rows); + RowVectorType rowvec = RowVectorType::Random(cols); + + // test addition + + m2 = m1; + m2.colwise() += colvec; + VERIFY_IS_APPROX(m2, m1.colwise() + colvec); + VERIFY_IS_APPROX(m2.col(c), m1.col(c) + colvec); + + VERIFY_RAISES_ASSERT(m2.colwise() += colvec.transpose()); + VERIFY_RAISES_ASSERT(m1.colwise() + colvec.transpose()); + + m2 = m1; + m2.rowwise() += rowvec; + VERIFY_IS_APPROX(m2, m1.rowwise() + rowvec); + VERIFY_IS_APPROX(m2.row(r), m1.row(r) + rowvec); + + VERIFY_RAISES_ASSERT(m2.rowwise() += rowvec.transpose()); + VERIFY_RAISES_ASSERT(m1.rowwise() + rowvec.transpose()); + + // test substraction + + m2 = m1; + m2.colwise() -= colvec; + VERIFY_IS_APPROX(m2, m1.colwise() - colvec); + VERIFY_IS_APPROX(m2.col(c), m1.col(c) - colvec); + + VERIFY_RAISES_ASSERT(m2.colwise() -= colvec.transpose()); + VERIFY_RAISES_ASSERT(m1.colwise() - colvec.transpose()); + + m2 = m1; + m2.rowwise() -= rowvec; + VERIFY_IS_APPROX(m2, m1.rowwise() - rowvec); + VERIFY_IS_APPROX(m2.row(r), m1.row(r) - rowvec); + + VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose()); + VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose()); +} + +void test_vectorwiseop() +{ + CALL_SUBTEST_1(vectorwiseop_array(Array22cd())); + CALL_SUBTEST_2(vectorwiseop_array(Array())); + CALL_SUBTEST_3(vectorwiseop_array(ArrayXXf(3, 4))); + CALL_SUBTEST_4(vectorwiseop_matrix(Matrix4cf())); + CALL_SUBTEST_5(vectorwiseop_matrix(Matrix())); + CALL_SUBTEST_6(vectorwiseop_matrix(MatrixXd(7,2))); +}