merge with eigen-tip
diff --git a/.hgignore b/.hgignore index 412e037..9432702 100644 --- a/.hgignore +++ b/.hgignore
@@ -20,4 +20,5 @@ tags .*.swp activity.png -gmon.out +*.out +*.php* \ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt index c229a28..0c068e6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt
@@ -33,20 +33,11 @@ set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) -option(EIGEN_BUILD_TESTS "Build tests" OFF) -option(EIGEN_BUILD_DEMOS "Build demos" OFF) -if(NOT WIN32) - option(EIGEN_BUILD_LIB "Build the binary shared library" OFF) -endif(NOT WIN32) option(EIGEN_BUILD_BTL "Build benchmark suite" OFF) if(NOT WIN32) option(EIGEN_BUILD_PKGCONFIG "Build pkg-config .pc file for Eigen" ON) endif(NOT WIN32) -if(EIGEN_BUILD_LIB) - option(EIGEN_TEST_LIB "Build the unit tests using the library (disable -pedantic)" OFF) -endif(EIGEN_BUILD_LIB) - set(CMAKE_INCLUDE_CURRENT_DIR ON) if(CMAKE_COMPILER_IS_GNUCXX) @@ -95,7 +86,7 @@ option(EIGEN_TEST_SSE2 "Enable/Disable SSE2 in tests/examples" OFF) if(EIGEN_TEST_SSE2) if(NOT CMAKE_CL_64) - # arch is not supported on 64 bit systems, SSE is enabled automatically. + # arch is not supported on 64 bit systems, SSE is enabled automatically. set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /arch:SSE2") endif(NOT CMAKE_CL_64) message("Enabling SSE2 in tests/examples") @@ -108,6 +99,10 @@ message("Disabling vectorization in tests/examples") endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION) +option(EIGEN_TEST_C++0x "Enables all C++0x features." OFF) + +option(EIGEN_TEST_RVALUE_REF_SUPPORT "Enable rvalue references for unit tests." OFF) + include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}) set(INCLUDE_INSTALL_DIR @@ -125,23 +120,43 @@ add_subdirectory(Eigen) -add_subdirectory(doc) +add_subdirectory(doc EXCLUDE_FROM_ALL) -if(EIGEN_BUILD_TESTS) - include(CTest) - add_subdirectory(test) -endif(EIGEN_BUILD_TESTS) +include(CTest) +enable_testing() # must be called from the root CMakeLists, see man page +if(EIGEN_CMAKE_RUN_FROM_CTEST) +add_subdirectory(test) # can't do EXCLUDE_FROM_ALL here, breaks CTest +else(EIGEN_CMAKE_RUN_FROM_CTEST) +add_subdirectory(test EXCLUDE_FROM_ALL) +endif(EIGEN_CMAKE_RUN_FROM_CTEST) add_subdirectory(unsupported) -if(EIGEN_BUILD_DEMOS) - add_subdirectory(demos) -endif(EIGEN_BUILD_DEMOS) +add_subdirectory(demos EXCLUDE_FROM_ALL) +add_subdirectory(blas EXCLUDE_FROM_ALL) + +# TODO: consider also replacing EIGEN_BUILD_BTL by a custom target "make btl"? if(EIGEN_BUILD_BTL) - add_subdirectory(bench/btl) + add_subdirectory(bench/btl EXCLUDE_FROM_ALL) endif(EIGEN_BUILD_BTL) -if(EIGEN_BUILD_TESTS) - ei_testing_print_summary() -endif(EIGEN_BUILD_TESTS) +ei_testing_print_summary() +if(NOT MSVC_IDE) +message("") +message("Configured Eigen ${EIGEN_VERSION_NUMBER}") +message("You can now do the following:") +message("--------------+----------------------------------------------------------------") +message("Command | Description") +message("--------------+----------------------------------------------------------------") +message("make install | Install to ${CMAKE_INSTALL_PREFIX}") +message(" | * To change that: cmake . -DCMAKE_INSTALL_PREFIX=yourpath") +message("make btest | Build the unit tests") +message(" | * That takes lots of memory! Easy on the -j option") +message("make test | Build and run the unit tests (using CTest)") +message("make test_qr | Build a specific test, here test_qr. To run it: test/test_qr") +message("make debug_qr | Build a test with full debug info. To run it: test/debug_qr") +message("make blas | Build BLAS library (not the same thing as Eigen)") +message("make doc | Generate the API documentation, requires Doxygen & LaTeX") +message("--------------+----------------------------------------------------------------") +endif(NOT MSVC_IDE)
diff --git a/Doxyfile b/Doxyfile index f31cf5b..6065beb 100644 --- a/Doxyfile +++ b/Doxyfile
@@ -5,7 +5,7 @@ #--------------------------------------------------------------------------- DOXYFILE_ENCODING = UTF-8 PROJECT_NAME = Eigen -PROJECT_NUMBER = 2.0 +PROJECT_NUMBER = you-got-it-wrong OUTPUT_DIRECTORY = ./ CREATE_SUBDIRS = NO OUTPUT_LANGUAGE = English
diff --git a/Eigen/Core b/Eigen/Core index 854f737..c8fcb1c 100644 --- a/Eigen/Core +++ b/Eigen/Core
@@ -106,6 +106,9 @@ #error The preprocessor symbols 'min' or 'max' are defined. If you are compiling on Windows, do #define NOMINMAX to prevent windows.h from defining these symbols. #endif +// defined in bits/termios.h +#undef B0 + namespace Eigen { /** \defgroup Core_Module Core module @@ -143,6 +146,7 @@ #include "src/Core/Functors.h" #include "src/Core/MatrixBase.h" +#include "src/Core/AnyMatrixBase.h" #include "src/Core/Coeffs.h" #ifndef EIGEN_PARSED_BY_DOXYGEN // work around Doxygen bug triggered by Assign.h r814874
diff --git a/Eigen/Householder b/Eigen/Householder index ba06bd8..ef3e613 100644 --- a/Eigen/Householder +++ b/Eigen/Householder
@@ -16,6 +16,7 @@ */ #include "src/Householder/Householder.h" +#include "src/Householder/HouseholderSequence.h" } // namespace Eigen
diff --git a/Eigen/src/Array/Replicate.h b/Eigen/src/Array/Replicate.h index 02f9c06..b20bcd4 100644 --- a/Eigen/src/Array/Replicate.h +++ b/Eigen/src/Array/Replicate.h
@@ -69,15 +69,22 @@ EIGEN_GENERIC_PUBLIC_INTERFACE(Replicate) - inline Replicate(const MatrixType& matrix) + template<typename OriginalMatrixType> + inline Replicate(const OriginalMatrixType& matrix) : m_matrix(matrix), m_rowFactor(RowFactor), m_colFactor(ColFactor) { + EIGEN_STATIC_ASSERT((ei_is_same_type<MatrixType,OriginalMatrixType>::ret), + THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE) ei_assert(RowFactor!=Dynamic && ColFactor!=Dynamic); } - inline Replicate(const MatrixType& matrix, int rowFactor, int colFactor) + template<typename OriginalMatrixType> + inline Replicate(const OriginalMatrixType& matrix, int rowFactor, int colFactor) : m_matrix(matrix), m_rowFactor(rowFactor), m_colFactor(colFactor) - {} + { + EIGEN_STATIC_ASSERT((ei_is_same_type<MatrixType,OriginalMatrixType>::ret), + THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE) + } inline int rows() const { return m_matrix.rows() * m_rowFactor.value(); } inline int cols() const { return m_matrix.cols() * m_colFactor.value(); }
diff --git a/Eigen/src/Cholesky/CholeskyInstantiations.cpp b/Eigen/src/Cholesky/CholeskyInstantiations.cpp deleted file mode 100644 index 92902f1..0000000 --- a/Eigen/src/Cholesky/CholeskyInstantiations.cpp +++ /dev/null
@@ -1,35 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> -// -// 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 <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_EXTERN_INSTANTIATIONS -#define EIGEN_EXTERN_INSTANTIATIONS -#endif -#include "../../Core" -#undef EIGEN_EXTERN_INSTANTIATIONS - -#include "../../Cholesky" - -namespace Eigen { - EIGEN_CHOLESKY_MODULE_INSTANTIATE(); -}
diff --git a/Eigen/src/Core/AnyMatrixBase.h b/Eigen/src/Core/AnyMatrixBase.h new file mode 100644 index 0000000..58b4257 --- /dev/null +++ b/Eigen/src/Core/AnyMatrixBase.h
@@ -0,0 +1,153 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_ANYMATRIXBASE_H +#define EIGEN_ANYMATRIXBASE_H + + +/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T). + * + * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase. + * + * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc. + * + * Notice that this class is trivial, it is only used to disambiguate overloaded functions. + */ +template<typename Derived> struct AnyMatrixBase +{ + typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType; + + Derived& derived() { return *static_cast<Derived*>(this); } + const Derived& derived() const { return *static_cast<const Derived*>(this); } + + /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ + inline int rows() const { return derived().rows(); } + /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ + inline int cols() const { return derived().cols(); } + + /** \internal Don't use it, but do the equivalent: \code dst = *this; \endcode */ + template<typename Dest> inline void evalTo(Dest& dst) const + { derived().evalTo(dst); } + + /** \internal Don't use it, but do the equivalent: \code dst += *this; \endcode */ + template<typename Dest> inline void addToDense(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + typename Dest::PlainMatrixType res(rows(),cols()); + evalTo(res); + dst += res; + } + + /** \internal Don't use it, but do the equivalent: \code dst -= *this; \endcode */ + template<typename Dest> inline void subToDense(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + typename Dest::PlainMatrixType res(rows(),cols()); + evalTo(res); + dst -= res; + } + + /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheRight(*this); \endcode */ + template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + dst = dst * this->derived(); + } + + /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheLeft(*this); \endcode */ + template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + dst = this->derived() * dst; + } + +}; + +/*************************************************************************** +* Implementation of matrix base methods +***************************************************************************/ + +/** Copies the generic expression \a other into *this. \returns a reference to *this. + * The expression must provide a (templated) evalTo(Derived& dst) const function + * which does the actual job. In practice, this allows any user to write its own + * special matrix without having to modify MatrixBase */ +template<typename Derived> +template<typename OtherDerived> +Derived& MatrixBase<Derived>::operator=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().evalTo(derived()); + return derived(); +} + +template<typename Derived> +template<typename OtherDerived> +Derived& MatrixBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().addToDense(derived()); + return derived(); +} + +template<typename Derived> +template<typename OtherDerived> +Derived& MatrixBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().subToDense(derived()); + return derived(); +} + +/** replaces \c *this by \c *this * \a other. + * + * \returns a reference to \c *this + */ +template<typename Derived> +template<typename OtherDerived> +inline Derived& +MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheRight(derived()); + return derived(); +} + +/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=() */ +template<typename Derived> +template<typename OtherDerived> +inline void MatrixBase<Derived>::applyOnTheRight(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheRight(derived()); +} + +/** replaces \c *this by \c *this * \a other. */ +template<typename Derived> +template<typename OtherDerived> +inline void MatrixBase<Derived>::applyOnTheLeft(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheLeft(derived()); +} + +#endif // EIGEN_ANYMATRIXBASE_H
diff --git a/Eigen/src/Core/CMakeLists.txt b/Eigen/src/Core/CMakeLists.txt index a555be7..2346fc2 100644 --- a/Eigen/src/Core/CMakeLists.txt +++ b/Eigen/src/Core/CMakeLists.txt
@@ -5,13 +5,6 @@ DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core COMPONENT Devel ) -FILE(GLOB Eigen_Core_Product_SRCS "products/*.h") - -INSTALL(FILES - ${Eigen_Core_Product_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/products COMPONENT Devel - ) - - +ADD_SUBDIRECTORY(products) ADD_SUBDIRECTORY(util) ADD_SUBDIRECTORY(arch)
diff --git a/Eigen/src/Core/CommaInitializer.h b/Eigen/src/Core/CommaInitializer.h index e86f47a..328c558 100644 --- a/Eigen/src/Core/CommaInitializer.h +++ b/Eigen/src/Core/CommaInitializer.h
@@ -116,6 +116,9 @@ int m_row; // current row id int m_col; // current col id int m_currentBlockRows; // current block height + +private: + CommaInitializer& operator=(const CommaInitializer&); }; /** \anchor MatrixBaseCommaInitRef
diff --git a/Eigen/src/Core/CoreInstantiations.cpp b/Eigen/src/Core/CoreInstantiations.cpp deleted file mode 100644 index 3c021a8..0000000 --- a/Eigen/src/Core/CoreInstantiations.cpp +++ /dev/null
@@ -1,47 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> -// -// 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 <http://www.gnu.org/licenses/>. - -#ifdef EIGEN_EXTERN_INSTANTIATIONS -#undef EIGEN_EXTERN_INSTANTIATIONS -#endif - -#include "../../Core" - -namespace Eigen -{ - -#define EIGEN_INSTANTIATE_PRODUCT(TYPE) \ -template static void ei_cache_friendly_product<TYPE>( \ - int _rows, int _cols, int depth, \ - bool _lhsRowMajor, const TYPE* _lhs, int _lhsStride, \ - bool _rhsRowMajor, const TYPE* _rhs, int _rhsStride, \ - bool resRowMajor, TYPE* res, int resStride) - -EIGEN_INSTANTIATE_PRODUCT(float); -EIGEN_INSTANTIATE_PRODUCT(double); -EIGEN_INSTANTIATE_PRODUCT(int); -EIGEN_INSTANTIATE_PRODUCT(std::complex<float>); -EIGEN_INSTANTIATE_PRODUCT(std::complex<double>); - -}
diff --git a/Eigen/src/Core/Cwise.h b/Eigen/src/Core/Cwise.h index fbf206d..4b14332 100644 --- a/Eigen/src/Core/Cwise.h +++ b/Eigen/src/Core/Cwise.h
@@ -178,6 +178,9 @@ protected: ExpressionTypeNested m_matrix; + + private: + Cwise& operator=(const Cwise&); }; /** \returns a Cwise wrapper of *this providing additional coefficient-wise operations
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h index a746959..4665fe0 100644 --- a/Eigen/src/Core/DiagonalMatrix.h +++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -51,7 +51,7 @@ DenseMatrixType toDenseMatrix() const { return derived(); } template<typename DenseDerived> - void evalToDense(MatrixBase<DenseDerived> &other) const; + void evalTo(MatrixBase<DenseDerived> &other) const; template<typename DenseDerived> void addToDense(MatrixBase<DenseDerived> &other) const { other.diagonal() += diagonal(); } @@ -72,7 +72,7 @@ template<typename Derived> template<typename DenseDerived> -void DiagonalBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const +void DiagonalBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const { other.setZero(); other.diagonal() = diagonal();
diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 0c68d74..cbaeb83 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h
@@ -351,6 +351,8 @@ EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const { return ei_pmul(a, ei_pset1(m_other)); } const Scalar m_other; +private: + ei_scalar_multiple_op& operator=(const ei_scalar_multiple_op&); }; template<typename Scalar> struct ei_functor_traits<ei_scalar_multiple_op<Scalar> > @@ -378,6 +380,8 @@ EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const { return ei_pmul(a, ei_pset1(m_other)); } const Scalar m_other; +private: + ei_scalar_quotient1_impl& operator=(const ei_scalar_quotient1_impl&); }; template<typename Scalar> struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,true> > @@ -423,6 +427,8 @@ EIGEN_STRONG_INLINE const Scalar operator() (int, int = 0) const { return m_other; } EIGEN_STRONG_INLINE const PacketScalar packetOp() const { return ei_pset1(m_other); } const Scalar m_other; +private: + ei_scalar_constant_op& operator=(const ei_scalar_constant_op&); }; template<typename Scalar> struct ei_functor_traits<ei_scalar_constant_op<Scalar> >
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 40edf4a..05469b3 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h
@@ -317,4 +317,34 @@ return a <= b || ei_isApprox(a, b, prec); } +/************** +*** bool *** +**************/ + +template<> inline bool precision<bool>() { return 0; } +inline bool ei_real(bool x) { return x; } +inline bool& ei_real_ref(bool& x) { return x; } +inline bool ei_imag(bool) { return 0; } +inline bool ei_conj(bool x) { return x; } +inline bool ei_abs(bool x) { return x; } +inline bool ei_abs2(bool x) { return x; } +inline bool ei_sqrt(bool x) { return x; } + +template<> inline bool ei_random() +{ + return (ei_random<int>(0,1) == 1); +} +inline bool ei_isMuchSmallerThan(bool a, bool, bool = precision<bool>()) +{ + return !a; +} +inline bool ei_isApprox(bool a, bool b, bool = precision<bool>()) +{ + return a == b; +} +inline bool ei_isApproxOrLessThan(bool a, bool b, bool = precision<bool>()) +{ + return int(a) <= int(b); +} + #endif // EIGEN_MATHFUNCTIONS_H
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index c08f124..027e6bb 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h
@@ -399,8 +399,8 @@ return Base::lazyAssign(other.derived()); } - template<typename OtherDerived,typename OtherEvalType> - EIGEN_STRONG_INLINE Matrix& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func) + template<typename OtherDerived> + EIGEN_STRONG_INLINE Matrix& operator=(const ReturnByValue<OtherDerived>& func) { resize(func.rows(), func.cols()); return Base::operator=(func); @@ -504,8 +504,8 @@ _set_noalias(other); } /** Copy constructor with in-place evaluation */ - template<typename OtherDerived,typename OtherEvalType> - EIGEN_STRONG_INLINE Matrix(const ReturnByValue<OtherDerived,OtherEvalType>& other) + template<typename OtherDerived> + EIGEN_STRONG_INLINE Matrix(const ReturnByValue<OtherDerived>& other) { _check_template_params(); resize(other.rows(), other.cols()); @@ -538,7 +538,7 @@ * data pointers. */ template<typename OtherDerived> - void swap(const MatrixBase<OtherDerived>& other); + void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other); /** \name Map * These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects, @@ -707,6 +707,8 @@ { static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other) { + if (_this.rows() == other.rows() && _this.cols() == other.cols()) return; + // Note: Here is space for improvement. Basically, for conservativeResize(int,int), // neither RowsAtCompileTime or ColsAtCompileTime must be Dynamic. If only one of the // dimensions is dynamic, one could use either conservativeResize(int rows, NoChange_t) or @@ -728,6 +730,8 @@ { static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other) { + if (_this.rows() == other.rows() && _this.cols() == other.cols()) return; + // segment(...) will check whether Derived/OtherDerived are vectors! typename MatrixBase<Derived>::PlainMatrixType tmp(other); const int common_size = std::min<int>(_this.size(),tmp.size()); @@ -756,7 +760,7 @@ template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols> template<typename OtherDerived> -inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other) +inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other) { enum { SwapPointers = ei_is_same_type<Matrix, OtherDerived>::ret && Base::SizeAtCompileTime==Dynamic }; ei_matrix_swap_impl<Matrix, OtherDerived, bool(SwapPointers)>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 25a0545..387c113 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h
@@ -26,44 +26,6 @@ #ifndef EIGEN_MATRIXBASE_H #define EIGEN_MATRIXBASE_H - -/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T). - * - * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase. - * - * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc. - * - * Notice that this class is trivial, it is only used to disambiguate overloaded functions. - */ -template<typename Derived> struct AnyMatrixBase -{ - typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType; - - Derived& derived() { return *static_cast<Derived*>(this); } - const Derived& derived() const { return *static_cast<const Derived*>(this); } - /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ - inline int rows() const { return derived().rows(); } - /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ - inline int cols() const { return derived().cols(); } - - template<typename Dest> inline void evalTo(Dest& dst) const - { derived().evalTo(dst); } - - template<typename Dest> inline void addToDense(Dest& dst) const - { - typename Dest::PlainMatrixType res(rows(),cols()); - evalToDense(res); - dst += res; - } - - template<typename Dest> inline void subToDense(Dest& dst) const - { - typename Dest::PlainMatrixType res(rows(),cols()); - evalToDense(res); - dst -= res; - } -}; - /** \class MatrixBase * * \brief Base class for all matrices, vectors, and expressions @@ -96,7 +58,6 @@ #endif // not EIGEN_PARSED_BY_DOXYGEN { public: - #ifndef EIGEN_PARSED_BY_DOXYGEN using ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar, typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>::operator*; @@ -301,24 +262,17 @@ */ Derived& operator=(const MatrixBase& other); - /** Copies the generic expression \a other into *this. \returns a reference to *this. - * The expression must provide a (templated) evalToDense(Derived& dst) const function - * which does the actual job. In practice, this allows any user to write its own - * special matrix without having to modify MatrixBase */ template<typename OtherDerived> - Derived& operator=(const AnyMatrixBase<OtherDerived> &other) - { other.derived().evalToDense(derived()); return derived(); } + Derived& operator=(const AnyMatrixBase<OtherDerived> &other); template<typename OtherDerived> - Derived& operator+=(const AnyMatrixBase<OtherDerived> &other) - { other.derived().addToDense(derived()); return derived(); } + Derived& operator+=(const AnyMatrixBase<OtherDerived> &other); template<typename OtherDerived> - Derived& operator-=(const AnyMatrixBase<OtherDerived> &other) - { other.derived().subToDense(derived()); return derived(); } + Derived& operator-=(const AnyMatrixBase<OtherDerived> &other); - template<typename OtherDerived,typename OtherEvalType> - Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func); + template<typename OtherDerived> + Derived& operator=(const ReturnByValue<OtherDerived>& func); #ifndef EIGEN_PARSED_BY_DOXYGEN /** Copies \a other into *this without evaluating other. \returns a reference to *this. */ @@ -436,6 +390,12 @@ template<typename OtherDerived> Derived& operator*=(const AnyMatrixBase<OtherDerived>& other); + template<typename OtherDerived> + void applyOnTheLeft(const AnyMatrixBase<OtherDerived>& other); + + template<typename OtherDerived> + void applyOnTheRight(const AnyMatrixBase<OtherDerived>& other); + template<typename DiagonalDerived> const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight> operator*(const DiagonalBase<DiagonalDerived> &diagonal) const; @@ -632,7 +592,7 @@ { return typename ei_eval<Derived>::type(derived()); } template<typename OtherDerived> - void swap(const MatrixBase<OtherDerived>& other); + void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other); template<unsigned int Added> const Flagged<Derived, Added, 0> marked() const; @@ -675,8 +635,11 @@ typename ei_traits<Derived>::Scalar minCoeff() const; typename ei_traits<Derived>::Scalar maxCoeff() const; - typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col = 0) const; - typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col = 0) const; + typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col) const; + typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col) const; + + typename ei_traits<Derived>::Scalar minCoeff(int* index) const; + typename ei_traits<Derived>::Scalar maxCoeff(int* index) const; template<typename BinaryOp> typename ei_result_of<BinaryOp(typename ei_traits<Derived>::Scalar)>::type @@ -739,8 +702,10 @@ const LU<PlainMatrixType> lu() const; const PartialLU<PlainMatrixType> partialLu() const; const PlainMatrixType inverse() const; - void computeInverse(PlainMatrixType *result) const; - bool computeInverseWithCheck( PlainMatrixType *result ) const; + template<typename ResultType> + void computeInverse(ResultType *result) const; + template<typename ResultType> + bool computeInverseWithCheck(ResultType *result ) const; Scalar determinant() const; /////////// Cholesky module ///////////
diff --git a/Eigen/src/Core/MatrixStorage.h b/Eigen/src/Core/MatrixStorage.h index f67095d..73b17e6 100644 --- a/Eigen/src/Core/MatrixStorage.h +++ b/Eigen/src/Core/MatrixStorage.h
@@ -29,32 +29,38 @@ struct ei_constructor_without_unaligned_array_assert {}; /** \internal - * Static array automatically aligned if the total byte size is a multiple of 16 and the matrix options require auto alignment + * Static array. If the MatrixOptions require auto-alignment, the array will be automatically aligned: + * to 16 bytes boundary if the total size is a multiple of 16 bytes. */ template <typename T, int Size, int MatrixOptions, - bool Align = (!(MatrixOptions&DontAlign)) && (((Size*sizeof(T))&0xf)==0) -> struct ei_matrix_array -{ - EIGEN_ALIGN_128 T array[Size]; - - ei_matrix_array() - { - #ifndef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT - ei_assert((reinterpret_cast<size_t>(array) & 0xf) == 0 - && "this assertion is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html **** READ THIS WEB PAGE !!! ****"); - #endif - } - - ei_matrix_array(ei_constructor_without_unaligned_array_assert) {} -}; - -template <typename T, int Size, int MatrixOptions> struct ei_matrix_array<T,Size,MatrixOptions,false> + int Alignment = (MatrixOptions&DontAlign) ? 0 + : (((Size*sizeof(T))%16)==0) ? 16 + : 0 > +struct ei_matrix_array { T array[Size]; ei_matrix_array() {} ei_matrix_array(ei_constructor_without_unaligned_array_assert) {} }; +#ifdef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT + #define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask) +#else + #define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask) \ + ei_assert((reinterpret_cast<size_t>(array) & sizemask) == 0 \ + && "this assertion is explained here: " \ + "http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html" \ + " **** READ THIS WEB PAGE !!! ****"); +#endif + +template <typename T, int Size, int MatrixOptions> +struct ei_matrix_array<T, Size, MatrixOptions, 16> +{ + EIGEN_ALIGN16 T array[Size]; + ei_matrix_array() { EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(0xf) } + ei_matrix_array(ei_constructor_without_unaligned_array_assert) {} +}; + /** \internal * * \class ei_matrix_storage
diff --git a/Eigen/src/Core/NestByValue.h b/Eigen/src/Core/NestByValue.h index 85a6727..94a8f80 100644 --- a/Eigen/src/Core/NestByValue.h +++ b/Eigen/src/Core/NestByValue.h
@@ -102,6 +102,9 @@ protected: const ExpressionType m_expression; + + private: + NestByValue& operator=(const NestByValue&); }; /** \returns an expression of the temporary version of *this.
diff --git a/Eigen/src/Core/NoAlias.h b/Eigen/src/Core/NoAlias.h index 66d8d83..dc32c2f 100644 --- a/Eigen/src/Core/NoAlias.h +++ b/Eigen/src/Core/NoAlias.h
@@ -73,6 +73,9 @@ protected: ExpressionType& m_expression; + + private: + NoAlias& operator=(const NoAlias&); }; /** \returns a pseudo expression of \c *this with an operator= assuming
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index e7227d4..cc75165 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h
@@ -135,7 +135,7 @@ { typedef typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested; typedef typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested; - typedef GeneralProduct<Lhs, Rhs, UnrolledProduct> Type; + typedef GeneralProduct<LhsNested, RhsNested, UnrolledProduct> Type; }; @@ -211,6 +211,9 @@ { ei_outer_product_selector<(int(Dest::Flags)&RowMajorBit) ? RowMajor : ColMajor>::run(*this, dest, alpha); } + + private: + GeneralProduct& operator=(const GeneralProduct&); }; template<> struct ei_outer_product_selector<ColMajor> { @@ -276,6 +279,9 @@ ei_gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor, bool(ei_blas_traits<MatrixType>::ActualAccess)>::run(*this, dst, alpha); } + +private: + GeneralProduct& operator=(const GeneralProduct&); }; // The vector is on the left => transposition @@ -434,18 +440,4 @@ return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); } - - -/** replaces \c *this by \c *this * \a other. - * - * \returns a reference to \c *this - */ -template<typename Derived> -template<typename OtherDerived> -inline Derived & -MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other) -{ - return derived() = derived() * other.derived(); -} - #endif // EIGEN_PRODUCT_H
diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index 764dc4d..57183ac 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h
@@ -137,6 +137,8 @@ void coeffRef(int,int); void coeff(int) const; void coeffRef(int); + + ProductBase& operator=(const ProductBase&); }; template<typename NestedProduct>
diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h index 3f2b478..4a5d5c1 100644 --- a/Eigen/src/Core/ReturnByValue.h +++ b/Eigen/src/Core/ReturnByValue.h
@@ -28,45 +28,37 @@ /** \class ReturnByValue * */ -template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int _MaxRows,int _MaxCols> -struct ei_traits<ReturnByValue<Functor,Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > > - : public ei_traits<Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > +template<typename Derived> +struct ei_traits<ReturnByValue<Derived> > + : public ei_traits<typename ei_traits<Derived>::ReturnMatrixType> { enum { - Flags = ei_traits<Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> >::Flags | EvalBeforeNestingBit + Flags = ei_traits<typename ei_traits<Derived>::ReturnMatrixType>::Flags | EvalBeforeNestingBit }; }; -template<typename Functor,typename EvalTypeDerived,int n> -struct ei_nested<ReturnByValue<Functor,MatrixBase<EvalTypeDerived> >, n, EvalTypeDerived> +template<typename Derived,int n,typename PlainMatrixType> +struct ei_nested<ReturnByValue<Derived>, n, PlainMatrixType> { - typedef EvalTypeDerived type; + typedef typename ei_traits<Derived>::ReturnMatrixType type; }; -template<typename Functor, typename EvalType> class ReturnByValue +template<typename Derived> + class ReturnByValue : public MatrixBase<ReturnByValue<Derived> > { - public: - template<typename Dest> inline void evalTo(Dest& dst) const - { static_cast<const Functor*>(this)->evalTo(dst); } -}; - -template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int _MaxRows,int _MaxCols> - class ReturnByValue<Functor,Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > - : public MatrixBase<ReturnByValue<Functor,Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > > -{ - typedef Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> EvalType; + typedef typename ei_traits<Derived>::ReturnMatrixType ReturnMatrixType; public: EIGEN_GENERIC_PUBLIC_INTERFACE(ReturnByValue) template<typename Dest> inline void evalTo(Dest& dst) const - { static_cast<const Functor* const>(this)->evalTo(dst); } - inline int rows() const { return static_cast<const Functor* const>(this)->rows(); } - inline int cols() const { return static_cast<const Functor* const>(this)->cols(); } + { static_cast<const Derived* const>(this)->evalTo(dst); } + inline int rows() const { return static_cast<const Derived* const>(this)->rows(); } + inline int cols() const { return static_cast<const Derived* const>(this)->cols(); } }; template<typename Derived> -template<typename OtherDerived,typename OtherEvalType> -Derived& MatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived,OtherEvalType>& other) +template<typename OtherDerived> +Derived& MatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived>& other) { other.evalTo(derived()); return derived();
diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 77fe797..06e69c4 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h
@@ -56,7 +56,7 @@ { const int blockSize = 4096; RealScalar scale = 0; - RealScalar invScale; + RealScalar invScale = 1; RealScalar ssq = 0; // sum of square enum { Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? ForceAligned : AsRequested @@ -108,21 +108,15 @@ iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent rbig = std::numeric_limits<RealScalar>::max(); // largest floating-point number - // Check the basic machine-dependent constants. - if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) - || (it<=4 && ibeta <= 3 ) || it<2) - { - ei_assert(false && "the algorithm cannot be guaranteed on this computer"); - } iexp = -((1-iemin)/2); - b1 = RealScalar(std::pow(double(ibeta),iexp)); // lower boundary of midrange + b1 = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // lower boundary of midrange iexp = (iemax + 1 - it)/2; - b2 = RealScalar(std::pow(double(ibeta),iexp)); // upper boundary of midrange + b2 = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // upper boundary of midrange iexp = (2-iemin)/2; - s1m = RealScalar(std::pow(double(ibeta),iexp)); // scaling factor for lower range + s1m = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for lower range iexp = - ((iemax+it)/2); - s2m = RealScalar(std::pow(double(ibeta),iexp)); // scaling factor for upper range + s2m = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range overfl = rbig*s2m; // overfow boundary for abig eps = RealScalar(std::pow(double(ibeta), 1-it));
diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h index 44e1f07..a7cf219 100644 --- a/Eigen/src/Core/Swap.h +++ b/Eigen/src/Core/Swap.h
@@ -128,15 +128,9 @@ */ template<typename Derived> template<typename OtherDerived> -void MatrixBase<Derived>::swap(const MatrixBase<OtherDerived>& other) +void MatrixBase<Derived>::swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other) { (SwapWrapper<Derived>(derived())).lazyAssign(other); } #endif // EIGEN_SWAP_H - - - - - -
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index b0362f2..e60d57e 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h
@@ -91,9 +91,9 @@ #endif // not EIGEN_PARSED_BY_DOXYGEN template<typename DenseDerived> - void evalToDense(MatrixBase<DenseDerived> &other) const; + void evalTo(MatrixBase<DenseDerived> &other) const; template<typename DenseDerived> - void evalToDenseLazy(MatrixBase<DenseDerived> &other) const; + void evalToLazy(MatrixBase<DenseDerived> &other) const; protected: @@ -300,13 +300,13 @@ } template<typename OtherDerived> - void swap(const TriangularBase<OtherDerived>& other) + void swap(TriangularBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other) { TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); } template<typename OtherDerived> - void swap(const MatrixBase<OtherDerived>& other) + void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other) { TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); } @@ -546,23 +546,23 @@ * If the matrix is triangular, the opposite part is set to zero. */ template<typename Derived> template<typename DenseDerived> -void TriangularBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const +void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const { if(ei_traits<Derived>::Flags & EvalBeforeAssigningBit) { typename Derived::PlainMatrixType other_evaluated(rows(), cols()); - evalToDenseLazy(other_evaluated); + evalToLazy(other_evaluated); other.derived().swap(other_evaluated); } else - evalToDenseLazy(other.derived()); + evalToLazy(other.derived()); } /** Assigns a triangular or selfadjoint matrix to a dense matrix. * If the matrix is triangular, the opposite part is set to zero. */ template<typename Derived> template<typename DenseDerived> -void TriangularBase<Derived>::evalToDenseLazy(MatrixBase<DenseDerived> &other) const +void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const { const bool unroll = DenseDerived::SizeAtCompileTime * Derived::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT;
diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h index b291f7b..65268b6 100644 --- a/Eigen/src/Core/VectorBlock.h +++ b/Eigen/src/Core/VectorBlock.h
@@ -77,11 +77,12 @@ typedef Block<VectorType, ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size, - PacketAccess> Base; + PacketAccess> _Base; enum { IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1 }; public: + _EIGEN_GENERIC_PUBLIC_INTERFACE(VectorBlock, _Base) using Base::operator=; using Base::operator+=;
diff --git a/Eigen/src/Core/Visitor.h b/Eigen/src/Core/Visitor.h index 598c2db..590efc7 100644 --- a/Eigen/src/Core/Visitor.h +++ b/Eigen/src/Core/Visitor.h
@@ -164,7 +164,7 @@ /** \returns the minimum of all coefficients of *this * and puts in *row and *col its location. * - * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff() + * \sa MatrixBase::minCoeff(int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff() */ template<typename Derived> typename ei_traits<Derived>::Scalar @@ -177,6 +177,22 @@ return minVisitor.res; } +/** \returns the minimum of all coefficients of *this + * and puts in *index its location. + * + * \sa MatrixBase::minCoeff(int*,int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff() + */ +template<typename Derived> +typename ei_traits<Derived>::Scalar +MatrixBase<Derived>::minCoeff(int* index) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + ei_min_coeff_visitor<Scalar> minVisitor; + this->visit(minVisitor); + *index = (RowsAtCompileTime==1) ? minVisitor.col : minVisitor.row; + return minVisitor.res; +} + /** \returns the maximum of all coefficients of *this * and puts in *row and *col its location. * @@ -193,5 +209,20 @@ return maxVisitor.res; } +/** \returns the maximum of all coefficients of *this + * and puts in *index its location. + * + * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::minCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::maxCoeff() + */ +template<typename Derived> +typename ei_traits<Derived>::Scalar +MatrixBase<Derived>::maxCoeff(int* index) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + ei_max_coeff_visitor<Scalar> maxVisitor; + this->visit(maxVisitor); + *index = (RowsAtCompileTime==1) ? maxVisitor.col : maxVisitor.row; + return maxVisitor.res; +} #endif // EIGEN_VISITOR_H
diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index a9c1620..1526a4b 100644 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h
@@ -265,14 +265,14 @@ template<> inline float ei_pfirst(const v4f& a) { - float EIGEN_ALIGN_128 af[4]; + float EIGEN_ALIGN16 af[4]; vec_st(a, 0, af); return af[0]; } template<> inline int ei_pfirst(const v4i& a) { - int EIGEN_ALIGN_128 ai[4]; + int EIGEN_ALIGN16 ai[4]; vec_st(a, 0, ai); return ai[0]; } @@ -373,7 +373,7 @@ inline int ei_predux_mul(const v4i& a) { - EIGEN_ALIGN_128 int aux[4]; + EIGEN_ALIGN16 int aux[4]; ei_pstore(aux, a); return aux[0] * aux[1] * aux[2] * aux[3]; }
diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 52e666d..eb1c2d3 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h
@@ -76,16 +76,14 @@ #ifdef __GNUC__ // Sometimes GCC implements _mm_set1_p* using multiple moves, -// that is inefficient :( +// that is inefficient :( (e.g., see ei_gemm_pack_rhs) template<> EIGEN_STRONG_INLINE Packet4f ei_pset1<float>(const float& from) { Packet4f res = _mm_set_ss(from); - asm("shufps $0, %[x], %[x]" : [x] "+x" (res) : ); - return res; + return _mm_shuffle_ps(res,res,0); } template<> EIGEN_STRONG_INLINE Packet2d ei_pset1<double>(const double& from) { Packet2d res = _mm_set_sd(from); - asm("unpcklpd %[x], %[x]" : [x] "+x" (res) : ); - return res; + return _mm_unpacklo_pd(res,res); } #else template<> EIGEN_STRONG_INLINE Packet4f ei_pset1<float>(const float& from) { return _mm_set1_ps(from); } @@ -361,7 +359,7 @@ // after some experiments, it is seems this is the fastest way to implement it // for GCC (eg., reusing ei_pmul is very slow !) // TODO try to call _mm_mul_epu32 directly - EIGEN_ALIGN_128 int aux[4]; + EIGEN_ALIGN16 int aux[4]; ei_pstore(aux, a); return (aux[0] * aux[1]) * (aux[2] * aux[3]);; } @@ -380,7 +378,7 @@ { // after some experiments, it is seems this is the fastest way to implement it // for GCC (eg., it does not like using std::min after the ei_pstore !!) - EIGEN_ALIGN_128 int aux[4]; + EIGEN_ALIGN16 int aux[4]; ei_pstore(aux, a); register int aux0 = aux[0]<aux[1] ? aux[0] : aux[1]; register int aux2 = aux[2]<aux[3] ? aux[2] : aux[3]; @@ -401,7 +399,7 @@ { // after some experiments, it is seems this is the fastest way to implement it // for GCC (eg., it does not like using std::min after the ei_pstore !!) - EIGEN_ALIGN_128 int aux[4]; + EIGEN_ALIGN16 int aux[4]; ei_pstore(aux, a); register int aux0 = aux[0]>aux[1] ? aux[0] : aux[1]; register int aux2 = aux[2]>aux[3] ? aux[2] : aux[3];
diff --git a/Eigen/src/Core/products/CMakeLists.txt b/Eigen/src/Core/products/CMakeLists.txt new file mode 100644 index 0000000..21fc94a --- /dev/null +++ b/Eigen/src/Core/products/CMakeLists.txt
@@ -0,0 +1,6 @@ +FILE(GLOB Eigen_Core_Product_SRCS "*.h") + +INSTALL(FILES + ${Eigen_Core_Product_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/products COMPONENT Devel + )
diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index 5787503..a18e5ef 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h
@@ -57,8 +57,7 @@ if(ConjugateRhs) alpha = ei_conj(alpha); -// std::cerr << "prod " << size << " " << rhs.size() << "\n"; - + typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename ei_packet_traits<Scalar>::type Packet; const int PacketSize = sizeof(Packet)/sizeof(Scalar); @@ -69,9 +68,9 @@ const int PeelAlignedMask = PacketSize*peels-1; // How many coeffs of the result do we have to skip to be aligned. - // Here we assume data are at least aligned on the base scalar type that is mandatory anyway. - const int alignedStart = ei_alignmentOffset(res,size); - const int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0; + // Here we assume data are at least aligned on the base scalar type. + int alignedStart = ei_alignmentOffset(res,size); + int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0; const int peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart; const int alignmentStep = PacketSize>1 ? (PacketSize - lhsStride % PacketSize) & PacketAlignedMask : 0; @@ -84,12 +83,18 @@ // find how many columns do we have to skip to be aligned with the result (if possible) int skipColumns = 0; - if (PacketSize>1) + // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats) + if( (size_t(lhs)%sizeof(RealScalar)) || (size_t(res)%sizeof(RealScalar)) ) + { + alignedSize = 0; + alignedStart = 0; + } + else if (PacketSize>1) { ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize); while (skipColumns<PacketSize && - alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize)) + alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize)) ++skipColumns; if (skipColumns==PacketSize) { @@ -263,6 +268,7 @@ ei_conj_helper<ConjugateLhs,ConjugateRhs> cj; + typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename ei_packet_traits<Scalar>::type Packet; const int PacketSize = sizeof(Packet)/sizeof(Scalar); @@ -274,9 +280,10 @@ const int size = rhsSize; // How many coeffs of the result do we have to skip to be aligned. - // Here we assume data are at least aligned on the base scalar type that is mandatory anyway. - const int alignedStart = ei_alignmentOffset(rhs, size); - const int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0; + // Here we assume data are at least aligned on the base scalar type + // if that's not the case then vectorization is discarded, see below. + int alignedStart = ei_alignmentOffset(rhs, size); + int alignedSize = PacketSize>1 ? alignedStart + ((size-alignedStart) & ~PacketAlignedMask) : 0; const int peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart; const int alignmentStep = PacketSize>1 ? (PacketSize - lhsStride % PacketSize) & PacketAlignedMask : 0; @@ -289,7 +296,13 @@ // find how many rows do we have to skip to be aligned with rhs (if possible) int skipRows = 0; - if (PacketSize>1) + // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats) + if( (size_t(lhs)%sizeof(RealScalar)) || (size_t(rhs)%sizeof(RealScalar)) ) + { + alignedSize = 0; + alignedStart = 0; + } + else if (PacketSize>1) { ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
diff --git a/Eigen/src/Core/products/SelfadjointRank2Update.h b/Eigen/src/Core/products/SelfadjointRank2Update.h index 64fcf86..69cf189 100644 --- a/Eigen/src/Core/products/SelfadjointRank2Update.h +++ b/Eigen/src/Core/products/SelfadjointRank2Update.h
@@ -38,10 +38,8 @@ static void run(Scalar* mat, int stride, const UType& u, const VType& v, Scalar alpha) { const int size = u.size(); -// std::cerr << "lower \n" << u.transpose() << "\n" << v.transpose() << "\n\n"; for (int i=0; i<size; ++i) { -// std::cerr << Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) += (alpha * ei_conj(u.coeff(i))) * v.end(size-i) + (alpha * ei_conj(v.coeff(i))) * u.end(size-i);
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index c5f27d8..65e5ce6 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -64,7 +64,7 @@ template<typename ExpressionType> class Cwise; template<typename ExpressionType> class WithFormat; template<typename MatrixType> struct CommaInitializer; -template<typename Functor, typename EvalType> class ReturnByValue; +template<typename Derived> class ReturnByValue; template<typename _Scalar, int Rows=Dynamic, int Cols=Dynamic, int Supers=Dynamic, int Subs=Dynamic, int Options=0> class BandMatrix; @@ -123,6 +123,7 @@ template<typename MatrixType, unsigned int Options = 0> class JacobiSVD; template<typename MatrixType, int UpLo = LowerTriangular> class LLT; template<typename MatrixType> class LDLT; +template<typename VectorsType, typename CoeffsType> class HouseholderSequence; template<typename Scalar> class PlanarRotation; // Geometry module:
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index ec8337e..706b301 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h
@@ -29,8 +29,8 @@ #undef minor #define EIGEN_WORLD_VERSION 2 -#define EIGEN_MAJOR_VERSION 0 -#define EIGEN_MINOR_VERSION 52 +#define EIGEN_MAJOR_VERSION 90 +#define EIGEN_MINOR_VERSION 0 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ @@ -202,26 +202,34 @@ #define EIGEN_ASM_COMMENT(X) #endif -/* EIGEN_ALIGN_128 forces data to be 16-byte aligned, EVEN if vectorization (EIGEN_VECTORIZE) is disabled, +/* EIGEN_ALIGN_TO_BOUNDARY(n) forces data to be n-byte aligned. This is used to satisfy SIMD requirements. + * However, we do that EVEN if vectorization (EIGEN_VECTORIZE) is disabled, * so that vectorization doesn't affect binary compatibility. * * If we made alignment depend on whether or not EIGEN_VECTORIZE is defined, it would be impossible to link * vectorized and non-vectorized code. */ #if !EIGEN_ALIGN - #define EIGEN_ALIGN_128 + #define EIGEN_ALIGN_TO_BOUNDARY(n) #elif (defined __GNUC__) - #define EIGEN_ALIGN_128 __attribute__((aligned(16))) + #define EIGEN_ALIGN_TO_BOUNDARY(n) __attribute__((aligned(n))) #elif (defined _MSC_VER) - #define EIGEN_ALIGN_128 __declspec(align(16)) + #define EIGEN_ALIGN_TO_BOUNDARY(n) __declspec(align(n)) #elif (defined __SUNPRO_CC) // FIXME not sure about this one: - #define EIGEN_ALIGN_128 __attribute__((aligned(16))) + #define EIGEN_ALIGN_TO_BOUNDARY(n) __attribute__((aligned(n))) #else - #error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler + #error Please tell me what is the equivalent of __attribute__((aligned(n))) for your compiler #endif -#define EIGEN_RESTRICT __restrict +#define EIGEN_ALIGN16 EIGEN_ALIGN_TO_BOUNDARY(16) + +#ifdef EIGEN_DONT_USE_RESTRICT_KEYWORD + #define EIGEN_RESTRICT +#endif +#ifndef EIGEN_RESTRICT + #define EIGEN_RESTRICT __restrict +#endif #ifndef EIGEN_STACK_ALLOCATION_LIMIT #define EIGEN_STACK_ALLOCATION_LIMIT 1000000 @@ -246,6 +254,13 @@ // needed to define it here as escaping characters in CMake add_definition's argument seems very problematic. #define EIGEN_DOCS_IO_FORMAT IOFormat(3, 0, " ", "\n", "", "") +// C++0x features +#if defined(__GXX_EXPERIMENTAL_CXX0X__) || (defined(_MSC_VER) && (_MSC_VER >= 1600)) + #define EIGEN_REF_TO_TEMPORARY && +#else + #define EIGEN_REF_TO_TEMPORARY const & +#endif + #ifdef _MSC_VER #define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \ using Base::operator =; \
diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index 73d302f..883f2d9 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h
@@ -41,7 +41,7 @@ #ifndef EIGEN_NO_STATIC_ASSERT - #ifdef __GXX_EXPERIMENTAL_CXX0X__ + #if defined(__GXX_EXPERIMENTAL_CXX0X__) || (defined(_MSC_VER) && (_MSC_VER >= 1600)) // if native static_assert is enabled, let's use it #define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG); @@ -77,7 +77,8 @@ THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES, INVALID_MATRIX_TEMPLATE_PARAMETERS, BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER, - THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX + THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX, + THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE }; };
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 6663819..86206ce 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -1,148 +1,149 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Claire Maurice -// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> -// -// 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 <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H -#define EIGEN_COMPLEX_EIGEN_SOLVER_H - -/** \eigenvalues_module \ingroup Eigenvalues_Module - * \nonstableyet - * - * \class ComplexEigenSolver - * - * \brief Eigen values/vectors solver for general complex matrices - * - * \param MatrixType the type of the matrix of which we are computing the eigen decomposition - * - * \sa class EigenSolver, class SelfAdjointEigenSolver - */ -template<typename _MatrixType> class ComplexEigenSolver -{ - public: - typedef _MatrixType MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; - typedef std::complex<RealScalar> Complex; - typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType; - typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType; - - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via ComplexEigenSolver::compute(const MatrixType&). - */ - ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) - {} - - ComplexEigenSolver(const MatrixType& matrix) - : m_eivec(matrix.rows(),matrix.cols()), - m_eivalues(matrix.cols()), - m_isInitialized(false) - { - compute(matrix); - } - - EigenvectorType eigenvectors(void) const - { - ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); - return m_eivec; - } - - EigenvalueType eigenvalues() const - { - ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); - return m_eivalues; - } - - void compute(const MatrixType& matrix); - - protected: - MatrixType m_eivec; - EigenvalueType m_eivalues; - bool m_isInitialized; -}; - - -template<typename MatrixType> -void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix) -{ - // this code is inspired from Jampack - assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - m_eivalues.resize(n,1); - - RealScalar eps = epsilon<RealScalar>(); - - // Reduce to complex Schur form - ComplexSchur<MatrixType> schur(matrix); - - m_eivalues = schur.matrixT().diagonal(); - - m_eivec.setZero(); - - Scalar d2, z; - RealScalar norm = matrix.norm(); - - // compute the (normalized) eigenvectors - for(int k=n-1 ; k>=0 ; k--) - { - d2 = schur.matrixT().coeff(k,k); - m_eivec.coeffRef(k,k) = Scalar(1.0,0.0); - for(int i=k-1 ; i>=0 ; i--) - { - m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k); - if(k-i-1>0) - m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value(); - z = schur.matrixT().coeff(i,i) - d2; - if(z==Scalar(0)) - ei_real_ref(z) = eps * norm; - m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z; - - } - m_eivec.col(k).normalize(); - } - - m_eivec = schur.matrixU() * m_eivec; - m_isInitialized = true; - - // sort the eigenvalues - { - for (int i=0; i<n; i++) - { - int k; - m_eivalues.cwise().abs().end(n-i).minCoeff(&k); - if (k != 0) - { - k += i; - std::swap(m_eivalues[k],m_eivalues[i]); - m_eivec.col(i).swap(m_eivec.col(k)); - } - } - } -} - - - -#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Claire Maurice +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H +#define EIGEN_COMPLEX_EIGEN_SOLVER_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class ComplexEigenSolver + * + * \brief Eigen values/vectors solver for general complex matrices + * + * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * + * \sa class EigenSolver, class SelfAdjointEigenSolver + */ +template<typename _MatrixType> class ComplexEigenSolver +{ + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> Complex; + typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType; + typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ComplexEigenSolver::compute(const MatrixType&). + */ + ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) + {} + + ComplexEigenSolver(const MatrixType& matrix) + : m_eivec(matrix.rows(),matrix.cols()), + m_eivalues(matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + EigenvectorType eigenvectors(void) const + { + ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); + return m_eivec; + } + + EigenvalueType eigenvalues() const + { + ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); + return m_eivalues; + } + + void compute(const MatrixType& matrix); + + protected: + MatrixType m_eivec; + EigenvalueType m_eivalues; + bool m_isInitialized; +}; + + +template<typename MatrixType> +void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix) +{ + // this code is inspired from Jampack + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + m_eivec.resize(n,n); + + RealScalar eps = epsilon<RealScalar>(); + + // Reduce to complex Schur form + ComplexSchur<MatrixType> schur(matrix); + + m_eivalues = schur.matrixT().diagonal(); + + m_eivec.setZero(); + + Scalar d2, z; + RealScalar norm = matrix.norm(); + + // compute the (normalized) eigenvectors + for(int k=n-1 ; k>=0 ; k--) + { + d2 = schur.matrixT().coeff(k,k); + m_eivec.coeffRef(k,k) = Scalar(1.0,0.0); + for(int i=k-1 ; i>=0 ; i--) + { + m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k); + if(k-i-1>0) + m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value(); + z = schur.matrixT().coeff(i,i) - d2; + if(z==Scalar(0)) + ei_real_ref(z) = eps * norm; + m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z; + + } + m_eivec.col(k).normalize(); + } + + m_eivec = schur.matrixU() * m_eivec; + m_isInitialized = true; + + // sort the eigenvalues + { + for (int i=0; i<n; i++) + { + int k; + m_eivalues.cwise().abs().end(n-i).minCoeff(&k); + if (k != 0) + { + k += i; + std::swap(m_eivalues[k],m_eivalues[i]); + m_eivec.col(i).swap(m_eivec.col(k)); + } + } + } +} + + + +#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 58e2ea4..0534715 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -31,8 +31,15 @@ * * \class ComplexShur * - * \brief Performs a complex Shur decomposition of a real or complex square matrix + * \brief Performs a complex Schur decomposition of a real or complex square matrix * + * Given a real or complex square matrix A, this class computes the Schur decomposition: + * \f$ A = U T U^*\f$ where U is a unitary complex matrix, and T is a complex upper + * triangular matrix. + * + * The diagonal of the matrix T corresponds to the eigenvalues of the matrix A. + * + * \sa class RealSchur, class EigenSolver */ template<typename _MatrixType> class ComplexSchur { @@ -42,41 +49,56 @@ typedef typename NumTraits<Scalar>::Real RealScalar; typedef std::complex<RealScalar> Complex; typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType; + enum { + Size = MatrixType::RowsAtCompileTime + }; - /** - * \brief Default Constructor. + /** \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to - * perform decompositions via ComplexSchur::compute(const MatrixType&). + * perform decompositions via ComplexSchur::compute(). */ - ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false) + ComplexSchur(int size = Size==Dynamic ? 0 : Size) + : m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false) {} - ComplexSchur(const MatrixType& matrix) + /** Constructor computing the Schur decomposition of the matrix \a matrix. + * If \a skipU is true, then the matrix U is not computed. */ + ComplexSchur(const MatrixType& matrix, bool skipU = false) : m_matT(matrix.rows(),matrix.cols()), m_matU(matrix.rows(),matrix.cols()), - m_isInitialized(false) + m_isInitialized(false), + m_matUisUptodate(false) { - compute(matrix); + compute(matrix, skipU); } - ComplexMatrixType matrixU() const + /** \returns a const reference to the matrix U of the respective Schur decomposition. */ + const ComplexMatrixType& matrixU() const { ei_assert(m_isInitialized && "ComplexSchur is not initialized."); + ei_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); return m_matU; } - ComplexMatrixType matrixT() const + /** \returns a const reference to the matrix T of the respective Schur decomposition. + * Note that this function returns a plain square matrix. If you want to reference + * only the upper triangular part, use: + * \code schur.matrixT().triangularView<Upper>() \endcode. */ + const ComplexMatrixType& matrixT() const { ei_assert(m_isInitialized && "ComplexShur is not initialized."); return m_matT; } - void compute(const MatrixType& matrix); + /** Computes the Schur decomposition of the matrix \a matrix. + * If \a skipU is true, then the matrix U is not computed. */ + void compute(const MatrixType& matrix, bool skipU = false); protected: ComplexMatrixType m_matT, m_matU; bool m_isInitialized; + bool m_matUisUptodate; }; /** Computes the principal value of the square root of the complex \a z. */ @@ -117,17 +139,20 @@ } template<typename MatrixType> -void ComplexSchur<MatrixType>::compute(const MatrixType& matrix) +void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) { // this code is inspired from Jampack + + m_matUisUptodate = false; assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); // Reduce to Hessenberg form + // TODO skip Q if skipU = true HessenbergDecomposition<MatrixType> hess(matrix); m_matT = hess.matrixH(); - m_matU = hess.matrixQ(); + if(!skipU) m_matU = hess.matrixQ(); int iu = m_matT.cols() - 1; int il; @@ -206,7 +231,7 @@ { m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint()); m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot); - m_matU.applyOnTheRight(i, i+1, rot); + if(!skipU) m_matU.applyOnTheRight(i, i+1, rot); if(i != iu-1) { @@ -232,6 +257,7 @@ */ m_isInitialized = true; + m_matUisUptodate = !skipU; } #endif // EIGEN_COMPLEX_SCHUR_H
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 3fc36c0..73d240d 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -194,6 +194,7 @@ assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); m_eivalues.resize(n,1); + m_eivec.resize(n,n); MatrixType matH = matrix; RealVectorType ort(n);
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index b1e21d4..bb7e3fc 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -88,14 +88,14 @@ _compute(m_matrix, m_hCoeffs); } - /** \returns the householder coefficients allowing to + /** \returns a const reference to the householder coefficients allowing to * reconstruct the matrix Q from the packed data. * * \sa packedMatrix() */ - CoeffVectorType householderCoefficients() const { return m_hCoeffs; } + const CoeffVectorType& householderCoefficients() const { return m_hCoeffs; } - /** \returns the internal result of the decomposition. + /** \returns a const reference to the internal representation of the decomposition. * * The returned matrix contains the following information: * - the upper part and lower sub-diagonal represent the Hessenberg matrix H
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 84856aa..9e155de 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -168,6 +168,7 @@ assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); m_eivalues.resize(n,1); + m_eivec.resize(n,n); if(n==1) {
diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index 2de99b5..035d213 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h
@@ -207,10 +207,19 @@ } template<typename MatrixType,typename Lhs> +struct ei_traits<ei_homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs> > +{ + typedef Matrix<typename ei_traits<MatrixType>::Scalar, + Lhs::RowsAtCompileTime, + MatrixType::ColsAtCompileTime, + MatrixType::PlainMatrixType::Options, + Lhs::MaxRowsAtCompileTime, + MatrixType::MaxColsAtCompileTime> ReturnMatrixType; +}; + +template<typename MatrixType,typename Lhs> struct ei_homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs> - : public ReturnByValue<ei_homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs>, - Matrix<typename ei_traits<MatrixType>::Scalar, - Lhs::RowsAtCompileTime,MatrixType::ColsAtCompileTime> > + : public ReturnByValue<ei_homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs> > { typedef typename ei_cleantype<typename Lhs::Nested>::type LhsNested; ei_homogeneous_left_product_impl(const Lhs& lhs, const MatrixType& rhs) @@ -236,10 +245,19 @@ }; template<typename MatrixType,typename Rhs> +struct ei_traits<ei_homogeneous_right_product_impl<Homogeneous<MatrixType,Horizontal>,Rhs> > +{ + typedef Matrix<typename ei_traits<MatrixType>::Scalar, + MatrixType::RowsAtCompileTime, + Rhs::ColsAtCompileTime, + MatrixType::PlainMatrixType::Options, + MatrixType::MaxRowsAtCompileTime, + Rhs::MaxColsAtCompileTime> ReturnMatrixType; +}; + +template<typename MatrixType,typename Rhs> struct ei_homogeneous_right_product_impl<Homogeneous<MatrixType,Horizontal>,Rhs> - : public ReturnByValue<ei_homogeneous_right_product_impl<Homogeneous<MatrixType,Horizontal>,Rhs>, - Matrix<typename ei_traits<MatrixType>::Scalar, - MatrixType::RowsAtCompileTime, Rhs::ColsAtCompileTime> > + : public ReturnByValue<ei_homogeneous_right_product_impl<Homogeneous<MatrixType,Horizontal>,Rhs> > { typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested; ei_homogeneous_right_product_impl(const MatrixType& lhs, const Rhs& rhs)
diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 9eb8ed5..d03fd52 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h
@@ -247,14 +247,14 @@ ei_transform_construct_from_matrix<OtherMatrixType,Mode,Dim,HDim>::run(this, other.matrix()); } - template<typename OtherDerived,typename OtherEvalType> - Transform(const ReturnByValue<OtherDerived,OtherEvalType>& other) + template<typename OtherDerived> + Transform(const ReturnByValue<OtherDerived>& other) { other.evalTo(*this); } - template<typename OtherDerived,typename OtherEvalType> - Transform& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& other) + template<typename OtherDerived> + Transform& operator=(const ReturnByValue<OtherDerived>& other) { other.evalTo(*this); return *this; @@ -345,6 +345,10 @@ /** \sa MatrixBase::setIdentity() */ void setIdentity() { m_matrix.setIdentity(); } + static const typename MatrixType::IdentityReturnType Identity() + { + return MatrixType::Identity(); + } template<typename OtherDerived> inline Transform& scale(const MatrixBase<OtherDerived> &other); @@ -395,7 +399,7 @@ Transform& fromPositionOrientationScale(const MatrixBase<PositionDerived> &position, const OrientationType& orientation, const MatrixBase<ScaleDerived> &scale); - inline const MatrixType inverse(TransformTraits traits = (TransformTraits)Mode) const; + inline Transform inverse(TransformTraits traits = (TransformTraits)Mode) const; /** \returns a const pointer to the column major internal matrix */ const Scalar* data() const { return m_matrix.data(); } @@ -874,7 +878,7 @@ /** \nonstableyet * - * \returns the inverse transformation matrix according to some given knowledge + * \returns the inverse transformation according to some given knowledge * on \c *this. * * \param traits allows to optimize the inversion process when the transformion @@ -892,37 +896,37 @@ * \sa MatrixBase::inverse() */ template<typename Scalar, int Dim, int Mode> -const typename Transform<Scalar,Dim,Mode>::MatrixType +Transform<Scalar,Dim,Mode> Transform<Scalar,Dim,Mode>::inverse(TransformTraits hint) const { + Transform res; if (hint == Projective) { - return m_matrix.inverse(); + res.matrix() = m_matrix.inverse(); } else { - MatrixType res; if (hint == Isometry) { - res.template corner<Dim,Dim>(TopLeft) = linear().transpose(); + res.matrix().template corner<Dim,Dim>(TopLeft) = linear().transpose(); } else if(hint&Affine) { - res.template corner<Dim,Dim>(TopLeft) = linear().inverse(); + res.matrix().template corner<Dim,Dim>(TopLeft) = linear().inverse(); } else { ei_assert(false && "Invalid transform traits in Transform::Inverse"); } // translation and remaining parts - res.template corner<Dim,1>(TopRight) = - res.template corner<Dim,Dim>(TopLeft) * translation(); + res.matrix().template corner<Dim,1>(TopRight) = - res.matrix().template corner<Dim,Dim>(TopLeft) * translation(); if(int(Mode)!=int(AffineCompact)) { - res.template block<1,Dim>(Dim,0).setZero(); - res.coeffRef(Dim,Dim) = 1; + res.matrix().template block<1,Dim>(Dim,0).setZero(); + res.matrix().coeffRef(Dim,Dim) = 1; } - return res; } + return res; } /*****************************************************
diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 36f02d7..775b2f8 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h
@@ -55,12 +55,12 @@ * \f$ H = I - tau v v^*\f$ * and the vector v is: * \f$ v^T = [1 essential^T] \f$ - * + * * On output: * \param essential the essential part of the vector \c v * \param tau the scaling factor of the householder transformation * \param beta the result of H * \c *this - * + * * \sa MatrixBase::makeHouseholderInPlace(), MatrixBase::applyHouseholderOnTheLeft(), * MatrixBase::applyHouseholderOnTheRight() */ @@ -73,10 +73,10 @@ { EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); - + RealScalar tailSqNorm = size()==1 ? 0 : tail.squaredNorm(); Scalar c0 = coeff(0); - + if(tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0)) { *tau = 0;
diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h new file mode 100644 index 0000000..f5a8dd7 --- /dev/null +++ b/Eigen/src/Householder/HouseholderSequence.h
@@ -0,0 +1,170 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H +#define EIGEN_HOUSEHOLDER_SEQUENCE_H + +/** \ingroup Householder_Module + * \householder_module + * \class HouseholderSequence + * \brief Represents a sequence of householder reflections with decreasing size + * + * This class represents a product sequence of householder reflections \f$ H = \Pi_0^{n-1} H_i \f$ + * where \f$ H_i \f$ is the i-th householder transformation \f$ I - h_i v_i v_i^* \f$, + * \f$ v_i \f$ is the i-th householder vector \f$ [ 1, m_vectors(i+1,i), m_vectors(i+2,i), ...] \f$ + * and \f$ h_i \f$ is the i-th householder coefficient \c m_coeffs[i]. + * + * Typical usages are listed below, where H is a HouseholderSequence: + * \code + * A.applyOnTheRight(H); // A = A * H + * A.applyOnTheLeft(H); // A = H * A + * A.applyOnTheRight(H.adjoint()); // A = A * H^* + * A.applyOnTheLeft(H.adjoint()); // A = H^* * A + * MatrixXd Q = H; // conversion to a dense matrix + * \endcode + * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate. + * + * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() + */ + +template<typename VectorsType, typename CoeffsType> +struct ei_traits<HouseholderSequence<VectorsType,CoeffsType> > +{ + typedef typename VectorsType::Scalar Scalar; + enum { + RowsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime, + ColsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime, + MaxRowsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime, + MaxColsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime, + Flags = 0 + }; +}; + +template<typename VectorsType, typename CoeffsType> class HouseholderSequence + : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> > +{ + typedef typename VectorsType::Scalar Scalar; + public: + + typedef HouseholderSequence<VectorsType, + typename ei_meta_if<NumTraits<Scalar>::IsComplex, + NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >, + CoeffsType>::ret> ConjugateReturnType; + + HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false) + : m_vectors(v), m_coeffs(h), m_trans(trans) + {} + + int rows() const { return m_vectors.rows(); } + int cols() const { return m_vectors.rows(); } + + HouseholderSequence transpose() const + { return HouseholderSequence(m_vectors, m_coeffs, !m_trans); } + + ConjugateReturnType conjugate() const + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); } + + ConjugateReturnType adjoint() const + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); } + + ConjugateReturnType inverse() const { return adjoint(); } + + /** \internal */ + template<typename DestType> void evalTo(DestType& dst) const + { + int vecs = std::min(m_vectors.cols(),m_vectors.rows()); + int length = m_vectors.rows(); + dst.setIdentity(); + Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows()); + for(int k = vecs-1; k >= 0; --k) + { + if(m_trans) + dst.corner(BottomRight, length-k, length-k) + .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0)); + else + dst.corner(BottomRight, length-k, length-k) + .applyHouseholderOnTheLeft(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(k)); + } + } + + /** \internal */ + template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const + { + int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector + Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows()); + for(int k = 0; k < vecs; ++k) + { + int actual_k = m_trans ? vecs-k-1 : k; + dst.corner(BottomRight, dst.rows(), length-k) + .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0)); + } + } + + /** \internal */ + template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const + { + int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector + Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols()); + for(int k = 0; k < vecs; ++k) + { + int actual_k = m_trans ? k : vecs-k-1; + dst.corner(BottomRight, length-actual_k, dst.cols()) + .applyHouseholderOnTheLeft(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); + } + } + + template<typename OtherDerived> + typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other) const + { + typename OtherDerived::PlainMatrixType res(other); + applyThisOnTheLeft(res); + return res; + } + + template<typename OtherDerived> friend + typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence& h) + { + typename OtherDerived::PlainMatrixType res(other); + h.applyThisOnTheRight(res); + return res; + } + + protected: + typename VectorsType::Nested m_vectors; + typename CoeffsType::Nested m_coeffs; + bool m_trans; + +private: + HouseholderSequence& operator=(const HouseholderSequence&); +}; + +template<typename VectorsType, typename CoeffsType> +HouseholderSequence<VectorsType,CoeffsType> makeHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false) +{ + return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans); +} + +#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 248b480..3255c10 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h
@@ -95,8 +95,8 @@ return true; } -template<typename MatrixType> -bool ei_compute_inverse_size4_helper(const MatrixType& matrix, MatrixType* result) +template<typename MatrixType, typename ResultType> +bool ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType* result) { /* Let's split M into four 2x2 blocks: * (P Q) @@ -195,47 +195,47 @@ *** Part 2 : selector and MatrixBase methods *** ***********************************************/ -template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime> +template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime> struct ei_compute_inverse { - static inline void run(const MatrixType& matrix, MatrixType* result) + static inline void run(const MatrixType& matrix, ResultType* result) { matrix.partialLu().computeInverse(result); } }; -template<typename MatrixType> -struct ei_compute_inverse<MatrixType, 1> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse<MatrixType, ResultType, 1> { - static inline void run(const MatrixType& matrix, MatrixType* result) + static inline void run(const MatrixType& matrix, ResultType* result) { typedef typename MatrixType::Scalar Scalar; result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0); } }; -template<typename MatrixType> -struct ei_compute_inverse<MatrixType, 2> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse<MatrixType, ResultType, 2> { - static inline void run(const MatrixType& matrix, MatrixType* result) + static inline void run(const MatrixType& matrix, ResultType* result) { ei_compute_inverse_size2(matrix, result); } }; -template<typename MatrixType> -struct ei_compute_inverse<MatrixType, 3> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse<MatrixType, ResultType, 3> { - static inline void run(const MatrixType& matrix, MatrixType* result) + static inline void run(const MatrixType& matrix, ResultType* result) { - ei_compute_inverse_size3<false, MatrixType, MatrixType>(matrix, result); + ei_compute_inverse_size3<false, MatrixType, ResultType>(matrix, result); } }; -template<typename MatrixType> -struct ei_compute_inverse<MatrixType, 4> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse<MatrixType, ResultType, 4> { - static inline void run(const MatrixType& matrix, MatrixType* result) + static inline void run(const MatrixType& matrix, ResultType* result) { ei_compute_inverse_size4_with_check(matrix, result); } @@ -256,11 +256,12 @@ * \sa inverse(), computeInverseWithCheck() */ template<typename Derived> -inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const +template<typename ResultType> +inline void MatrixBase<Derived>::computeInverse(ResultType *result) const { ei_assert(rows() == cols()); EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) - ei_compute_inverse<PlainMatrixType>::run(eval(), result); + ei_compute_inverse<PlainMatrixType, ResultType>::run(eval(), result); } /** \lu_module @@ -291,10 +292,10 @@ * Compute inverse with invertibility check * *******************************************/ -template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime> +template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime> struct ei_compute_inverse_with_check { - static inline bool run(const MatrixType& matrix, MatrixType* result) + static inline bool run(const MatrixType& matrix, ResultType* result) { typedef typename MatrixType::Scalar Scalar; LU<MatrixType> lu( matrix ); @@ -304,10 +305,10 @@ } }; -template<typename MatrixType> -struct ei_compute_inverse_with_check<MatrixType, 1> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse_with_check<MatrixType, ResultType, 1> { - static inline bool run(const MatrixType& matrix, MatrixType* result) + static inline bool run(const MatrixType& matrix, ResultType* result) { typedef typename MatrixType::Scalar Scalar; if( matrix.coeff(0,0) == Scalar(0) ) return false; @@ -316,28 +317,28 @@ } }; -template<typename MatrixType> -struct ei_compute_inverse_with_check<MatrixType, 2> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse_with_check<MatrixType, ResultType, 2> { - static inline bool run(const MatrixType& matrix, MatrixType* result) + static inline bool run(const MatrixType& matrix, ResultType* result) { return ei_compute_inverse_size2_with_check(matrix, result); } }; -template<typename MatrixType> -struct ei_compute_inverse_with_check<MatrixType, 3> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse_with_check<MatrixType, ResultType, 3> { - static inline bool run(const MatrixType& matrix, MatrixType* result) + static inline bool run(const MatrixType& matrix, ResultType* result) { - return ei_compute_inverse_size3<true, MatrixType, MatrixType>(matrix, result); + return ei_compute_inverse_size3<true, MatrixType, ResultType>(matrix, result); } }; -template<typename MatrixType> -struct ei_compute_inverse_with_check<MatrixType, 4> +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse_with_check<MatrixType, ResultType, 4> { - static inline bool run(const MatrixType& matrix, MatrixType* result) + static inline bool run(const MatrixType& matrix, ResultType* result) { return ei_compute_inverse_size4_with_check(matrix, result); } @@ -354,11 +355,12 @@ * \sa inverse(), computeInverse() */ template<typename Derived> -inline bool MatrixBase<Derived>::computeInverseWithCheck(PlainMatrixType *result) const +template<typename ResultType> +inline bool MatrixBase<Derived>::computeInverseWithCheck(ResultType *result) const { ei_assert(rows() == cols()); EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) - return ei_compute_inverse_with_check<PlainMatrixType>::run(eval(), result); + return ei_compute_inverse_with_check<PlainMatrixType, ResultType>::run(eval(), result); }
diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h index 0ef59ba..e467c62 100644 --- a/Eigen/src/LU/PartialLU.h +++ b/Eigen/src/LU/PartialLU.h
@@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -215,10 +216,10 @@ typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU; typedef Block<MapLU, Dynamic, Dynamic> MatrixType; typedef Block<MatrixType,Dynamic,Dynamic> BlockType; - + /** \internal performs the LU decomposition in-place of the matrix \a lu * using an unblocked algorithm. - * + * * In addition, this function returns the row transpositions in the * vector \a row_transpositions which must have a size equal to the number * of columns of the matrix \a lu, and an integer \a nb_transpositions @@ -232,7 +233,7 @@ for(int k = 0; k < size; ++k) { int row_of_biggest_in_col; - lu.block(k,k,rows-k,1).cwise().abs().maxCoeff(&row_of_biggest_in_col); + lu.col(k).end(rows-k).cwise().abs().maxCoeff(&row_of_biggest_in_col); row_of_biggest_in_col += k; row_transpositions[k] = row_of_biggest_in_col; @@ -295,7 +296,7 @@ int bs = std::min(size-k,blockSize); // actual size of the block int trows = rows - k - bs; // trailing rows int tsize = size - k - bs; // trailing size - + // partition the matrix: // A00 | A01 | A02 // lu = A10 | A11 | A12 @@ -343,7 +344,7 @@ { ei_assert(lu.cols() == row_transpositions.size()); ei_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); - + ei_partial_lu_impl <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.stride(), &row_transpositions.coeffRef(0), nb_transpositions);
diff --git a/Eigen/src/QR/ColPivotingHouseholderQR.h b/Eigen/src/QR/ColPivotingHouseholderQR.h index 8024e3b..b141da0 100644 --- a/Eigen/src/QR/ColPivotingHouseholderQR.h +++ b/Eigen/src/QR/ColPivotingHouseholderQR.h
@@ -45,14 +45,14 @@ template<typename MatrixType> class ColPivotingHouseholderQR { public: - + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; @@ -62,6 +62,7 @@ typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType; + typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; /** * \brief Default Constructor. @@ -99,7 +100,7 @@ template<typename OtherDerived, typename ResultType> bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; - MatrixQType matrixQ(void) const; + HouseholderSequenceType matrixQ(void) const; /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -110,13 +111,13 @@ } ColPivotingHouseholderQR& compute(const MatrixType& matrix); - + const IntRowVectorType& colsPermutation() const { ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); return m_cols_permutation; } - + /** \returns the absolute value of the determinant of the matrix of which * *this is the QR decomposition. It has only linear complexity * (that is, O(n) where n is the dimension of the square matrix) @@ -145,7 +146,7 @@ * \sa absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; - + /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This is computed at the time of the construction of the QR decomposition. This @@ -268,7 +269,7 @@ int cols = matrix.cols(); int size = std::min(rows,cols); m_rank = size; - + m_qr = matrix; m_hCoeffs.resize(size); @@ -279,18 +280,18 @@ IntRowVectorType cols_transpositions(matrix.cols()); m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - + RealRowVectorType colSqNorms(cols); for(int k = 0; k < cols; ++k) colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); RealScalar biggestColSqNorm = colSqNorms.maxCoeff(); - + for (int k = 0; k < size; ++k) { int biggest_col_in_corner; RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner); biggest_col_in_corner += k; - + // if the corner is negligible, then we have less than full rank, and we can finish early if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision)) { @@ -302,10 +303,11 @@ } break; } - + cols_transpositions.coeffRef(k) = biggest_col_in_corner; if(k != biggest_col_in_corner) { m_qr.col(k).swap(m_qr.col(biggest_col_in_corner)); + std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_in_corner)); ++number_of_transpositions; } @@ -315,7 +317,7 @@ m_qr.corner(BottomRight, rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); - + colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2(); } @@ -325,7 +327,7 @@ m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -349,18 +351,12 @@ } const int rows = m_qr.rows(); - const int cols = b.cols(); ei_assert(b.rows() == rows); - + typename OtherDerived::PlainMatrixType c(b); - - Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols); - for (int k = 0; k < m_rank; ++k) - { - int remainingSize = rows-k; - c.corner(BottomRight, remainingSize, cols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); - } + + // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T + c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose()); if(!isSurjective()) { @@ -380,25 +376,12 @@ return true; } -/** \returns the matrix Q */ +/** \returns the matrix Q as a sequence of householder transformations */ template<typename MatrixType> -typename ColPivotingHouseholderQR<MatrixType>::MatrixQType ColPivotingHouseholderQR<MatrixType>::matrixQ() const +typename ColPivotingHouseholderQR<MatrixType>::HouseholderSequenceType ColPivotingHouseholderQR<MatrixType>::matrixQ() const { ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); - // compute the product H'_0 H'_1 ... H'_n-1, - // where H_k is the k-th Householder transformation I - h_k v_k v_k' - // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] - int rows = m_qr.rows(); - int cols = m_qr.cols(); - int size = std::min(rows,cols); - MatrixQType res = MatrixQType::Identity(rows, rows); - Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows); - for (int k = size-1; k >= 0; k--) - { - res.block(k, k, rows-k, rows-k) - .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); - } - return res; + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); } #endif // EIGEN_HIDE_HEAVY_CODE
diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index 0d542cf..9fee778 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h
@@ -45,14 +45,14 @@ template<typename MatrixType> class FullPivotingHouseholderQR { public: - + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; @@ -106,13 +106,13 @@ } FullPivotingHouseholderQR& compute(const MatrixType& matrix); - + const IntRowVectorType& colsPermutation() const { ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); return m_cols_permutation; } - + const IntColVectorType& rowsTranspositions() const { ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); @@ -147,7 +147,7 @@ * \sa absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; - + /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This is computed at the time of the construction of the QR decomposition. This @@ -271,7 +271,7 @@ int cols = matrix.cols(); int size = std::min(rows,cols); m_rank = size; - + m_qr = matrix; m_hCoeffs.resize(size); @@ -283,9 +283,9 @@ IntRowVectorType cols_transpositions(matrix.cols()); m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - + RealScalar biggest(0); - + for (int k = 0; k < size; ++k) { int row_of_biggest_in_corner, col_of_biggest_in_corner; @@ -297,7 +297,7 @@ row_of_biggest_in_corner += k; col_of_biggest_in_corner += k; if(k==0) biggest = biggest_in_corner; - + // if the corner is negligible, then we have less than full rank, and we can finish early if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) { @@ -336,7 +336,7 @@ m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -358,13 +358,13 @@ } else return false; } - + const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); - + typename OtherDerived::PlainMatrixType c(b); - + Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols); for (int k = 0; k < m_rank; ++k) {
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index a893058..4cc5539 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h
@@ -56,12 +56,13 @@ Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; + typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType; typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; + typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; /** * \brief Default Constructor. @@ -97,7 +98,12 @@ template<typename OtherDerived, typename ResultType> void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; - MatrixQType matrixQ(void) const; + MatrixQType matrixQ() const; + + HouseholderSequenceType matrixQAsHouseholderSequence() const + { + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); + } /** \returns a reference to the matrix where the Householder QR decomposition is stored * in a LAPACK-compatible way. @@ -169,7 +175,7 @@ int rows = matrix.rows(); int cols = matrix.cols(); int size = std::min(rows,cols); - + m_qr = matrix; m_hCoeffs.resize(size); @@ -200,26 +206,22 @@ ) const { ei_assert(m_isInitialized && "HouseholderQR is not initialized."); + result->resize(m_qr.cols(), b.cols()); const int rows = m_qr.rows(); - const int cols = b.cols(); + const int rank = std::min(m_qr.rows(), m_qr.cols()); ei_assert(b.rows() == rows); - result->resize(rows, cols); - *result = b; - - Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols); - for (int k = 0; k < cols; ++k) - { - int remainingSize = rows-k; + typename OtherDerived::PlainMatrixType c(b); - result->corner(BottomRight, remainingSize, cols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); - } + // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T + c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,rank), m_hCoeffs.start(rank)).transpose()); - const int rank = std::min(result->rows(), result->cols()); m_qr.corner(TopLeft, rank, rank) .template triangularView<UpperTriangular>() - .solveInPlace(result->corner(TopLeft, rank, result->cols())); + .solveInPlace(c.corner(TopLeft, rank, c.cols())); + + result->corner(TopLeft, rank, c.cols()) = c.corner(TopLeft,rank, c.cols()); + result->corner(BottomLeft, result->rows()-rank, c.cols()).setZero(); } /** \returns the matrix Q */ @@ -227,20 +229,7 @@ typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const { ei_assert(m_isInitialized && "HouseholderQR is not initialized."); - // compute the product H'_0 H'_1 ... H'_n-1, - // where H_k is the k-th Householder transformation I - h_k v_k v_k' - // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] - int rows = m_qr.rows(); - int cols = m_qr.cols(); - int size = std::min(rows,cols); - MatrixQType res = MatrixQType::Identity(rows, rows); - Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows); - for (int k = size-1; k >= 0; k--) - { - res.block(k, k, rows-k, rows-k) - .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); - } - return res; + return matrixQAsHouseholderSequence(); } #endif // EIGEN_HIDE_HEAVY_CODE
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 2801ee0..6a05978 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h
@@ -25,6 +25,23 @@ #ifndef EIGEN_JACOBISVD_H #define EIGEN_JACOBISVD_H +// forward declarations (needed by ICC) +// the empty bodies are required by VC +template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> +struct ei_svd_precondition_2x2_block_to_be_real {}; + +template<typename MatrixType, unsigned int Options, + bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0 + && (MatrixType::RowsAtCompileTime==Dynamic + || (MatrixType::RowsAtCompileTime>MatrixType::ColsAtCompileTime))> +struct ei_svd_precondition_if_more_rows_than_cols; + +template<typename MatrixType, unsigned int Options, + bool PossiblyMoreColsThanRows = (Options & AtLeastAsManyRowsAsCols) == 0 + && (MatrixType::ColsAtCompileTime==Dynamic + || (MatrixType::ColsAtCompileTime>MatrixType::RowsAtCompileTime))> +struct ei_svd_precondition_if_more_cols_than_rows; + /** \ingroup SVD_Module * \nonstableyet * @@ -118,8 +135,8 @@ friend struct ei_svd_precondition_if_more_cols_than_rows; }; -template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> -struct ei_svd_precondition_2x2_block_to_be_real +template<typename MatrixType, unsigned int Options> +struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options, false> { typedef JacobiSVD<MatrixType, Options> SVD; static void run(typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&, int, int) {} @@ -195,10 +212,7 @@ *j_left = rot1 * j_right->transpose(); } -template<typename MatrixType, unsigned int Options, - bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0 - && (MatrixType::RowsAtCompileTime==Dynamic - || MatrixType::RowsAtCompileTime>MatrixType::ColsAtCompileTime)> +template<typename MatrixType, unsigned int Options, bool PossiblyMoreRowsThanCols> struct ei_svd_precondition_if_more_rows_than_cols { typedef JacobiSVD<MatrixType, Options> SVD; @@ -231,10 +245,7 @@ } }; -template<typename MatrixType, unsigned int Options, - bool PossiblyMoreColsThanRows = (Options & AtLeastAsManyRowsAsCols) == 0 - && (MatrixType::ColsAtCompileTime==Dynamic - || MatrixType::ColsAtCompileTime>MatrixType::RowsAtCompileTime)> +template<typename MatrixType, unsigned int Options, bool PossiblyMoreColsThanRows> struct ei_svd_precondition_if_more_cols_than_rows { typedef JacobiSVD<MatrixType, Options> SVD; @@ -256,7 +267,7 @@ MaxColsAtCompileTime = SVD::MaxColsAtCompileTime, MatrixOptions = SVD::MatrixOptions }; - + static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd) { int rows = matrix.rows();
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index d0b7c98..61e8ade 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -452,7 +452,7 @@ /** \internal use operator= */ template<typename DenseDerived> - void evalToDense(MatrixBase<DenseDerived>& dst) const + void evalTo(MatrixBase<DenseDerived>& dst) const { dst.setZero(); for (int j=0; j<outerSize(); ++j) @@ -515,7 +515,7 @@ { return typename ei_eval<Derived>::type(derived()); } // template<typename OtherDerived> -// void swap(const MatrixBase<OtherDerived>& other); +// void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other); template<unsigned int Added> const SparseFlagged<Derived, Added, 0> marked() const;
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h index a9ad8aa..cb1196b 100644 --- a/Eigen/src/Sparse/SparseProduct.h +++ b/Eigen/src/Sparse/SparseProduct.h
@@ -305,23 +305,9 @@ } // dense = sparse * dense -// template<typename Derived> -// template<typename Lhs, typename Rhs> -// Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product) -// { -// typedef typename ei_cleantype<Lhs>::type _Lhs; -// typedef typename _Lhs::InnerIterator LhsInnerIterator; -// enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit }; -// derived().setZero(); -// for (int j=0; j<product.lhs().outerSize(); ++j) -// for (LhsInnerIterator i(product.lhs(),j); i; ++i) -// derived().row(LhsIsRowMajor ? j : i.index()) += i.value() * product.rhs().row(LhsIsRowMajor ? i.index() : j); -// return derived(); -// } - -template<typename Derived> -template<typename Lhs, typename Rhs> -Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product) +// Note that here we force no inlining and separate the setZero() because GCC messes up otherwise +template<typename Lhs, typename Rhs, typename Dest> +EIGEN_DONT_INLINE void ei_sparse_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) { typedef typename ei_cleantype<Lhs>::type _Lhs; typedef typename ei_cleantype<Rhs>::type _Rhs; @@ -335,34 +321,44 @@ || ( (_Lhs::Flags&LowerTriangularBit) && LhsIsRowMajor) ), ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf) }; - derived().setZero(); - for (int j=0; j<product.lhs().outerSize(); ++j) + for (int j=0; j<lhs.outerSize(); ++j) { - LhsInnerIterator i(product.lhs(),j); + LhsInnerIterator i(lhs,j); if (ProcessSecondHalf && i && (i.index()==j)) { - derived().row(j) += i.value() * product.rhs().row(j); + dst.row(j) += i.value() * rhs.row(j); ++i; } - Block<Derived,1,Derived::ColsAtCompileTime> res(derived().row(LhsIsRowMajor ? j : 0)); - for (; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) + typename Rhs::Scalar rhs_j = rhs.coeff(j,0); + Block<Dest,1,Dest::ColsAtCompileTime> dst_j(dst.row(LhsIsRowMajor ? j : 0)); + for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) { - if (LhsIsSelfAdjoint) + if(LhsIsSelfAdjoint) { int a = LhsIsRowMajor ? j : i.index(); int b = LhsIsRowMajor ? i.index() : j; - Scalar v = i.value(); - derived().row(a) += (v) * product.rhs().row(b); - derived().row(b) += ei_conj(v) * product.rhs().row(a); + typename Lhs::Scalar v = i.value(); + dst.row(a) += (v) * rhs.row(b); + dst.row(b) += ei_conj(v) * rhs.row(a); } - else if (LhsIsRowMajor) - res += i.value() * product.rhs().row(i.index()); + else if(LhsIsRowMajor) + dst_j += i.value() * rhs.row(i.index()); + else if(Rhs::ColsAtCompileTime==1) + dst.coeffRef(i.index()) += i.value() * rhs_j; else - derived().row(i.index()) += i.value() * product.rhs().row(j); + dst.row(i.index()) += i.value() * rhs.row(j); } if (ProcessFirstHalf && i && (i.index()==j)) - derived().row(j) += i.value() * product.rhs().row(j); + dst.row(j) += i.value() * rhs.row(j); } +} + +template<typename Derived> +template<typename Lhs, typename Rhs> +Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product) +{ + derived().setZero(); + ei_sparse_time_dense_product(product.lhs(), product.rhs(), derived()); return derived(); }
diff --git a/INSTALL b/INSTALL new file mode 100644 index 0000000..4f717e9 --- /dev/null +++ b/INSTALL
@@ -0,0 +1,35 @@ +Installation instructions for Eigen +*********************************** + +Explanation before starting +*************************** + +Eigen consists only of header files, hence there is nothing to compile +before you can use it. Moreover, these header files do not depend on your +platform, they are the same for everybody. + +Method 1. Installing without using CMake +**************************************** + +You can use right away the headers in the Eigen/ subdirectory. In order +to install, just copy this Eigen/ subdirectory to your favorite location. +If you also want the unsupported features, copy the unsupported/ +subdirectory too. + +Method 2. Installing using CMake +******************************** + +Let's call this directory 'source_dir' (where this INSTALL file is). +Before starting, create another directory which we will call 'build_dir'. + +Do: + + cd build_dir + cmake source_dir + make install + +The "make install" step may require administrator privileges. + +You can adjust the installation destination (the "prefix") +by passing the -DCMAKE_INSTALL_PREFIX=myprefix option to cmake, as is +explained in the message that cmake prints at the end.
diff --git a/Mainpage.dox b/Mainpage.dox index b5b7f30..7e89501 100644 --- a/Mainpage.dox +++ b/Mainpage.dox
@@ -1,127 +1,19 @@ - -// Please don't remove the following lines: -// this is the only way to specify doxygen options -// to api.kde.org's scripts - -// DOXYGEN_SET_PROJECT_NAME = Eigen2 -// DOXYGEN_SET_PROJECT_NUMBER = "2.0 - trunk" - -// DOXYGEN_SET_CREATE_SUBDIRS = NO -// DOXYGEN_SET_BRIEF_MEMBER_DESC = YES -// DOXYGEN_SET_REPEAT_BRIEF = YES -// DOXYGEN_SET_ALWAYS_DETAILED_SEC = NO -// DOXYGEN_SET_INLINE_INHERITED_MEMB = NO -// DOXYGEN_SET_FULL_PATH_NAMES = NO -// DOXYGEN_SET_SHORT_NAMES = NO -// DOXYGEN_SET_JAVADOC_AUTOBRIEF = NO -// DOXYGEN_SET_QT_AUTOBRIEF = NO -// DOXYGEN_SET_MULTILINE_CPP_IS_BRIEF = NO -// DOXYGEN_SET_DETAILS_AT_TOP = YES -// DOXYGEN_SET_INHERIT_DOCS = YES -// DOXYGEN_SET_ALIASES = "only_for_vectors=This is only for vectors (either row-vectors or column-vectors), i.e. matrices which are known at compile-time to have either one row or one column." "array_module=This is defined in the %Array module. \code #include <Eigen/Array> \endcode" "lu_module=This is defined in the %LU module. \code #include <Eigen/LU> \endcode" "cholesky_module=This is defined in the %Cholesky module. \code #include <Eigen/Cholesky> \endcode" "qr_module=This is defined in the %QR module. \code #include <Eigen/QR> \endcode" "svd_module=This is defined in the %SVD module. \code #include <Eigen/SVD> \endcode" "geometry_module=This is defined in the %Geometry module. \code #include <Eigen/Geometry> \endcode" "leastsquares_module=This is defined in the %LeastSquares module. \code #include <Eigen/LeastSquares> \endcode" "addexample=\anchor" "label=\bug" "redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" -// DOXYGEN_SET_DISTRIBUTE_GROUP_DOC = NO -// DOXYGEN_SET_SUBGROUPING = YES -// DOXYGEN_SET_TYPEDEF_HIDES_STRUCT = NO - -// DOXYGEN_SET_EXTRACT_ALL = NO -// DOXYGEN_SET_EXTRACT_PRIVATE = NO -// DOXYGEN_SET_EXTRACT_STATIC = NO -// DOXYGEN_SET_EXTRACT_LOCAL_CLASSES = NO -// DOXYGEN_SET_EXTRACT_LOCAL_METHODS = NO -// DOXYGEN_SET_EXTRACT_ANON_NSPACES = NO -// DOXYGEN_SET_HIDE_UNDOC_MEMBERS = NO -// DOXYGEN_SET_HIDE_UNDOC_CLASSES = YES -// DOXYGEN_SET_HIDE_FRIEND_COMPOUNDS = YES -// DOXYGEN_SET_HIDE_IN_BODY_DOCS = NO -// DOXYGEN_SET_INTERNAL_DOCS = NO -// DOXYGEN_SET_CASE_SENSE_NAMES = YES -// DOXYGEN_SET_HIDE_SCOPE_NAMES = YES -// DOXYGEN_SET_SHOW_INCLUDE_FILES = YES -// DOXYGEN_SET_INLINE_INFO = YES -// DOXYGEN_SET_SORT_MEMBER_DOCS = YES -// DOXYGEN_SET_SORT_BRIEF_DOCS = YES -// DOXYGEN_SET_SORT_GROUP_NAMES = NO -// DOXYGEN_SET_SORT_BY_SCOPE_NAME = NO -// DOXYGEN_SET_GENERATE_TODOLIST = NO -// DOXYGEN_SET_GENERATE_TESTLIST = NO -// DOXYGEN_SET_GENERATE_BUGLIST = NO -// DOXYGEN_SET_GENERATE_DEPRECATEDLIST= NO -// DOXYGEN_SET_SHOW_USED_FILES = YES -// DOXYGEN_SET_SHOW_DIRECTORIES = NO -// DOXYGEN_SET_SHOW_FILES = YES -// DOXYGEN_SET_SHOW_NAMESPACES = NO - -// DOXYGEN_SET_WARN_IF_UNDOCUMENTED = NO -// DOXYGEN_SET_WARN_NO_PARAMDOC = NO - -// DOXYGEN_SET_INPUT = @topdir@/eigen2/Eigen @topdir@/eigen2/doc @topdir@/eigen2/build/doc @topdir@/eigen2/unsupported/Eigen -// DOXYGEN_SET_EXCLUDE = *.sh *.in - -// DOXYGEN_SET_EXAMPLE_PATH = @topdir@/eigen2/doc/snippets/ @topdir@/eigen2/doc/examples/ @topdir@/eigen2/build/doc/examples/ @topdir@/eigen2/build/doc/snippets/ - -// DOXYGEN_SET_FILE_PATTERNS = * -// DOXYGEN_SET_RECURSIVE = NO -// DOXYGEN_SET_FILTER_SOURCE_FILES = YES - -// DOXYGEN_EXCLUDE_SYMBOLS = MatrixBase<* MapBase<* RotationBase<* Matrix<* - -// DOXYGEN_SET_SOURCE_BROWSER = NO -// DOXYGEN_SET_INLINE_SOURCES = NO -// DOXYGEN_SET_STRIP_CODE_COMMENTS = YES -// DOXYGEN_SET_REFERENCED_BY_RELATION = YES -// DOXYGEN_SET_REFERENCES_RELATION = YES -// DOXYGEN_SET_REFERENCES_LINK_SOURCE = YES -// DOXYGEN_SET_VERBATIM_HEADERS = YES - -// DOXYGEN_SET_ALPHABETICAL_INDEX = NO - -// DOXYGEN_SET_HTML_ALIGN_MEMBERS = YES -// DOXYGEN_SET_GENERATE_TREEVIEW = NO -// DOXYGEN_SET_FORMULA_FONTSIZE = 12 - -// DOXYGEN_SET_GENERATE_LATEX = NO -// DOXYGEN_SET_EXTRA_PACKAGES = amssymb - -// DOXYGEN_SET_ENABLE_PREPROCESSING = YES -// DOXYGEN_SET_MACRO_EXPANSION = YES -// DOXYGEN_SET_EXPAND_ONLY_PREDEF = YES -// DOXYGEN_SET_SEARCH_INCLUDES = YES -// DOXYGEN_SET_PREDEFINED = EIGEN_EMPTY_STRUCT EIGEN_PARSED_BY_DOXYGEN EIGEN_VECTORIZE EIGEN_QT_SUPPORT EIGEN_STRONG_INLINE=inline -// DOXYGEN_SET_EXPAND_AS_DEFINED = EIGEN_MAKE_SCALAR_OPS EIGEN_MAKE_TYPEDEFS EIGEN_MAKE_TYPEDEFS_ALL_SIZES EIGEN_CWISE_UNOP_RETURN_TYPE EIGEN_CWISE_BINOP_RETURN_TYPE -// DOXYGEN_SET_SKIP_FUNCTION_MACROS = YES - -// DOXYGEN_SET_CLASS_DIAGRAMS = NO -// DOXYGEN_SET_HIDE_UNDOC_RELATIONS = NO -// DOXYGEN_SET_HAVE_DOT = NO -// DOXYGEN_SET_CLASS_GRAPH = NO -// DOXYGEN_SET_COLLABORATION_GRAPH = NO -// DOXYGEN_SET_GROUP_GRAPHS = NO -// DOXYGEN_SET_UML_LOOK = NO -// DOXYGEN_SET_TEMPLATE_RELATIONS = NO -// DOXYGEN_SET_INCLUDE_GRAPH = NO -// DOXYGEN_SET_INCLUDED_BY_GRAPH = NO -// DOXYGEN_SET_CALL_GRAPH = NO -// DOXYGEN_SET_CALLER_GRAPH = NO -// DOXYGEN_SET_GRAPHICAL_HIERARCHY = NO -// DOXYGEN_SET_DIRECTORY_GRAPH = NO - - o /** \mainpage Eigen <h3>If you see this page, then you have not properly generated the documentation. Namely, you have run doxygen from the source directory, which is not appropriate for generating the documentation of Eigen.</h3> In order to generate the documentation of Eigen, please follow these steps: <ul> - <li>make sure you have the required software installed: cmake, doxygen, and a C++ compiler. + <li>make sure you have the required software installed: CMake, Doxygen, LaTeX, and a C++ compiler. <li>create a new directory, which we will call the "build directory", outside of the Eigen source directory.</li> <li>enter the build directory</li> <li>configure the project: <pre>cmake /path/to/source/directory</pre></li> - <li>now generate the documentaion: <pre>make doc</pre> or, if you have two CPUs, <pre>make doc -j2</pre> Note that this will compile the examples, run them, and integrate their output into the documentation, which can take some time.</li> + <li>now generate the documentaion: <pre>make doc</pre> or, if you have two CPUs, <pre>make doc -j3</pre> Note that this will compile the examples, run them, and integrate their output into the documentation, which can take some time.</li> </ul> After doing that, you will find the HTML documentation in the doc/html/ subdirectory of the build directory. -<h2>Note however that the documentation is available online here: -<a href="http://eigen.tuxfamily.org/dox/">http://eigen.tuxfamily.org/dox</a></h2> - +<h2>Note however that the documentation is available online here:</h2> +<a href="http://eigen.tuxfamily.org/dox/">http://eigen.tuxfamily.org/dox</a> for the stable version and +<a href="http://eigen.tuxfamily.org/dox-devel/">http://eigen.tuxfamily.org/dox-devel</a> for the development branch version. */
diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h index f4b67cb..f8dc8bd 100644 --- a/bench/BenchSparseUtil.h +++ b/bench/BenchSparseUtil.h
@@ -27,23 +27,23 @@ void fillMatrix(float density, int rows, int cols, EigenSparseMatrix& dst) { - dst.startFill(rows*cols*density); + dst.reserve(rows*cols*density); for(int j = 0; j < cols; j++) { for(int i = 0; i < rows; i++) { Scalar v = (ei_random<float>(0,1) < density) ? ei_random<Scalar>() : 0; if (v!=0) - dst.fillrand(i,j) = v; + dst.insert(i,j) = v; } } - dst.endFill(); + dst.finalize(); } void fillMatrix2(int nnzPerCol, int rows, int cols, EigenSparseMatrix& dst) { std::cout << "alloc " << nnzPerCol*cols << "\n"; - dst.startFill(nnzPerCol*cols); + dst.reserve(nnzPerCol*cols); for(int j = 0; j < cols; j++) { std::set<int> aux; @@ -54,10 +54,10 @@ k = ei_random<int>(0,rows-1); aux.insert(k); - dst.fillrand(k,j) = ei_random<Scalar>(); + dst.insert(k,j) = ei_random<Scalar>(); } } - dst.endFill(); + dst.finalize(); } void eiToDense(const EigenSparseMatrix& src, DenseMatrix& dst)
diff --git a/bench/sparse_dense_product.cpp b/bench/sparse_dense_product.cpp index 9b9af88..bfe4612 100644 --- a/bench/sparse_dense_product.cpp +++ b/bench/sparse_dense_product.cpp
@@ -91,22 +91,22 @@ { std::cout << "Eigen sparse\t" << sm1.nonZeros()/float(sm1.rows()*sm1.cols())*100 << "%\n"; - BENCH(for (int k=0; k<REPEAT; ++k) v2 = sm1 * v1;) + BENCH(asm("#myc"); v2 = sm1 * v1; asm("#myd");) std::cout << " a * v:\t" << timer.value() << endl; - - BENCH(for (int k=0; k<REPEAT; ++k) { asm("#mya"); v2 = sm1.transpose() * v1; asm("#myb"); }) - + + BENCH( { asm("#mya"); v2 = sm1.transpose() * v1; asm("#myb"); }) + std::cout << " a' * v:\t" << timer.value() << endl; } - + // { // DynamicSparseMatrix<Scalar> m1(sm1); // std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/float(m1.rows()*m1.cols())*100 << "%\n"; -// +// // BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1 * v1;) // std::cout << " a * v:\t" << timer.value() << endl; -// +// // BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1.transpose() * v1;) // std::cout << " a' * v:\t" << timer.value() << endl; // } @@ -118,23 +118,15 @@ //GmmDynSparse gmmT3(rows,cols); GmmSparse m1(rows,cols); eiToGmm(sm1, m1); - + std::vector<Scalar> gmmV1(cols), gmmV2(cols); Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1; Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2; - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - gmm::mult(m1, gmmV1, gmmV2); - timer.stop(); + BENCH( asm("#myx"); gmm::mult(m1, gmmV1, gmmV2); asm("#myy"); ) std::cout << " a * v:\t" << timer.value() << endl; - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - gmm::mult(gmm::transposed(m1), gmmV1, gmmV2); - timer.stop(); + BENCH( gmm::mult(gmm::transposed(m1), gmmV1, gmmV2); ) std::cout << " a' * v:\t" << timer.value() << endl; } #endif
diff --git a/bench/sparse_setter.cpp b/bench/sparse_setter.cpp index 6f7a19d..9c22636 100644 --- a/bench/sparse_setter.cpp +++ b/bench/sparse_setter.cpp
@@ -12,7 +12,15 @@ #endif #ifndef REPEAT -#define REPEAT 1 +#define REPEAT 2 +#endif + +#ifndef NBTRIES +#define NBTRIES 2 +#endif + +#ifndef KK +#define KK 10 #endif #ifndef NOGOOGLE @@ -22,7 +30,7 @@ #include "BenchSparseUtil.h" -#define CHECK_MEM +#define CHECK_MEM // #define CHECK_MEM std/**/::cout << "check mem\n"; getchar(); #define BENCH(X) \ @@ -37,9 +45,13 @@ typedef std::vector<float> Values; EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals); +EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals); +EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals); +EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals); EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals); EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals); EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals); +EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals); EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals); EIGEN_DONT_INLINE Scalar* setrand_ublas_coord(const Coordinates& coords, const Values& vals); EIGEN_DONT_INLINE Scalar* setrand_ublas_compressed(const Coordinates& coords, const Values& vals); @@ -50,17 +62,36 @@ { int rows = SIZE; int cols = SIZE; - bool fullyrand = false; - //float density = float(NBPERROW)/float(SIZE); - + bool fullyrand = true; + BenchTimer timer; Coordinates coords; Values values; if(fullyrand) { - for (int i=0; i<cols*NBPERROW; ++i) + Coordinates pool; + pool.reserve(cols*NBPERROW); + std::cerr << "fill pool" << "\n"; + for (int i=0; i<cols*NBPERROW; ) { - coords.push_back(Vector2i(ei_random<int>(0,rows-1),ei_random<int>(0,cols-1))); +// DynamicSparseMatrix<int> stencil(SIZE,SIZE); + Vector2i ij(ei_random<int>(0,rows-1),ei_random<int>(0,cols-1)); +// if(stencil.coeffRef(ij.x(), ij.y())==0) + { +// stencil.coeffRef(ij.x(), ij.y()) = 1; + pool.push_back(ij); + + } + ++i; + } + std::cerr << "pool ok" << "\n"; + int n = cols*NBPERROW*KK; + coords.reserve(n); + values.reserve(n); + for (int i=0; i<n; ++i) + { + int i = ei_random<int>(0,pool.size()); + coords.push_back(pool[i]); values.push_back(ei_random<Scalar>()); } } @@ -79,67 +110,55 @@ // dense matrices #ifdef DENSEMATRIX { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setrand_eigen_dense(coords,values); - timer.stop(); + BENCH(setrand_eigen_dense(coords,values);) std::cout << "Eigen Dense\t" << timer.value() << "\n"; } #endif // eigen sparse matrices - if (!fullyrand) +// if (!fullyrand) +// { +// BENCH(setinnerrand_eigen(coords,values);) +// std::cout << "Eigen fillrand\t" << timer.value() << "\n"; +// } { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setinnerrand_eigen(coords,values); - timer.stop(); - std::cout << "Eigen fillrand\t" << timer.value() << "\n"; + BENCH(setrand_eigen_dynamic(coords,values);) + std::cout << "Eigen dynamic\t" << timer.value() << "\n"; + } +// { +// BENCH(setrand_eigen_compact(coords,values);) +// std::cout << "Eigen compact\t" << timer.value() << "\n"; +// } + { + BENCH(setrand_eigen_sumeq(coords,values);) + std::cout << "Eigen sumeq\t" << timer.value() << "\n"; } { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setrand_eigen_gnu_hash(coords,values); - timer.stop(); - std::cout << "Eigen std::map\t" << timer.value() << "\n"; +// BENCH(setrand_eigen_gnu_hash(coords,values);) +// std::cout << "Eigen std::map\t" << timer.value() << "\n"; + } + { + BENCH(setrand_scipy(coords,values);) + std::cout << "scipy\t" << timer.value() << "\n"; } #ifndef NOGOOGLE { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setrand_eigen_google_dense(coords,values); - timer.stop(); + BENCH(setrand_eigen_google_dense(coords,values);) std::cout << "Eigen google dense\t" << timer.value() << "\n"; } { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setrand_eigen_google_sparse(coords,values); - timer.stop(); + BENCH(setrand_eigen_google_sparse(coords,values);) std::cout << "Eigen google sparse\t" << timer.value() << "\n"; } #endif - + #ifndef NOUBLAS { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setrand_ublas_mapped(coords,values); - timer.stop(); - std::cout << "ublas mapped\t" << timer.value() << "\n"; +// BENCH(setrand_ublas_mapped(coords,values);) +// std::cout << "ublas mapped\t" << timer.value() << "\n"; } { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setrand_ublas_genvec(coords,values); - timer.stop(); + BENCH(setrand_ublas_genvec(coords,values);) std::cout << "ublas vecofvec\t" << timer.value() << "\n"; } /*{ @@ -159,16 +178,12 @@ std::cout << "ublas coord\t" << timer.value() << "\n"; }*/ #endif - - + + // MTL4 #ifndef NOMTL { - timer.reset(); - timer.start(); - for (int k=0; k<REPEAT; ++k) - setrand_mtl(coords,values); - timer.stop(); + BENCH(setrand_mtl(coords,values)); std::cout << "MTL\t" << timer.value() << "\n"; } #endif @@ -180,16 +195,63 @@ { using namespace Eigen; SparseMatrix<Scalar> mat(SIZE,SIZE); - mat.startFill(2000000/*coords.size()*/); + //mat.startFill(2000000/*coords.size()*/); for (int i=0; i<coords.size(); ++i) { - mat.fillrand(coords[i].x(), coords[i].y()) = vals[i]; + mat.insert(coords[i].x(), coords[i].y()) = vals[i]; } - mat.endFill(); + mat.finalize(); CHECK_MEM; return 0; } +EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals) +{ + using namespace Eigen; + DynamicSparseMatrix<Scalar> mat(SIZE,SIZE); + mat.reserve(coords.size()/10); + for (int i=0; i<coords.size(); ++i) + { + mat.coeffRef(coords[i].x(), coords[i].y()) += vals[i]; + } + mat.finalize(); + CHECK_MEM; + return &mat.coeffRef(coords[0].x(), coords[0].y()); +} + +EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals) +{ + using namespace Eigen; + int n = coords.size()/KK; + DynamicSparseMatrix<Scalar> mat(SIZE,SIZE); + for (int j=0; j<KK; ++j) + { + DynamicSparseMatrix<Scalar> aux(SIZE,SIZE); + mat.reserve(n); + for (int i=j*n; i<(j+1)*n; ++i) + { + aux.insert(coords[i].x(), coords[i].y()) += vals[i]; + } + aux.finalize(); + mat += aux; + } + return &mat.coeffRef(coords[0].x(), coords[0].y()); +} + +EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals) +{ + using namespace Eigen; + DynamicSparseMatrix<Scalar> setter(SIZE,SIZE); + setter.reserve(coords.size()/10); + for (int i=0; i<coords.size(); ++i) + { + setter.coeffRef(coords[i].x(), coords[i].y()) += vals[i]; + } + SparseMatrix<Scalar> mat = setter; + CHECK_MEM; + return &mat.coeffRef(coords[0].x(), coords[0].y()); +} + EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals) { using namespace Eigen; @@ -198,11 +260,11 @@ RandomSetter<SparseMatrix<Scalar>, StdMapTraits > setter(mat); for (int i=0; i<coords.size(); ++i) { - setter(coords[i].x(), coords[i].y()) = vals[i]; + setter(coords[i].x(), coords[i].y()) += vals[i]; } CHECK_MEM; } - return 0;//&mat.coeffRef(coords[0].x(), coords[0].y()); + return &mat.coeffRef(coords[0].x(), coords[0].y()); } #ifndef NOGOOGLE @@ -213,10 +275,10 @@ { RandomSetter<SparseMatrix<Scalar>, GoogleDenseHashMapTraits> setter(mat); for (int i=0; i<coords.size(); ++i) - setter(coords[i].x(), coords[i].y()) = vals[i]; + setter(coords[i].x(), coords[i].y()) += vals[i]; CHECK_MEM; } - return 0;//&mat.coeffRef(coords[0].x(), coords[0].y()); + return &mat.coeffRef(coords[0].x(), coords[0].y()); } EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals) @@ -226,13 +288,139 @@ { RandomSetter<SparseMatrix<Scalar>, GoogleSparseHashMapTraits> setter(mat); for (int i=0; i<coords.size(); ++i) - setter(coords[i].x(), coords[i].y()) = vals[i]; + setter(coords[i].x(), coords[i].y()) += vals[i]; CHECK_MEM; } - return 0;//&mat.coeffRef(coords[0].x(), coords[0].y()); + return &mat.coeffRef(coords[0].x(), coords[0].y()); } #endif + +template <class T> +void coo_tocsr(const int n_row, + const int n_col, + const int nnz, + const Coordinates Aij, + const Values Ax, + int Bp[], + int Bj[], + T Bx[]) +{ + //compute number of non-zero entries per row of A coo_tocsr + std::fill(Bp, Bp + n_row, 0); + + for (int n = 0; n < nnz; n++){ + Bp[Aij[n].x()]++; + } + + //cumsum the nnz per row to get Bp[] + for(int i = 0, cumsum = 0; i < n_row; i++){ + int temp = Bp[i]; + Bp[i] = cumsum; + cumsum += temp; + } + Bp[n_row] = nnz; + + //write Aj,Ax into Bj,Bx + for(int n = 0; n < nnz; n++){ + int row = Aij[n].x(); + int dest = Bp[row]; + + Bj[dest] = Aij[n].y(); + Bx[dest] = Ax[n]; + + Bp[row]++; + } + + for(int i = 0, last = 0; i <= n_row; i++){ + int temp = Bp[i]; + Bp[i] = last; + last = temp; + } + + //now Bp,Bj,Bx form a CSR representation (with possible duplicates) +} + +template< class T1, class T2 > +bool kv_pair_less(const std::pair<T1,T2>& x, const std::pair<T1,T2>& y){ + return x.first < y.first; +} + + +template<class I, class T> +void csr_sort_indices(const I n_row, + const I Ap[], + I Aj[], + T Ax[]) +{ + std::vector< std::pair<I,T> > temp; + + for(I i = 0; i < n_row; i++){ + I row_start = Ap[i]; + I row_end = Ap[i+1]; + + temp.clear(); + + for(I jj = row_start; jj < row_end; jj++){ + temp.push_back(std::make_pair(Aj[jj],Ax[jj])); + } + + std::sort(temp.begin(),temp.end(),kv_pair_less<I,T>); + + for(I jj = row_start, n = 0; jj < row_end; jj++, n++){ + Aj[jj] = temp[n].first; + Ax[jj] = temp[n].second; + } + } +} + +template <class I, class T> +void csr_sum_duplicates(const I n_row, + const I n_col, + I Ap[], + I Aj[], + T Ax[]) +{ + I nnz = 0; + I row_end = 0; + for(I i = 0; i < n_row; i++){ + I jj = row_end; + row_end = Ap[i+1]; + while( jj < row_end ){ + I j = Aj[jj]; + T x = Ax[jj]; + jj++; + while( jj < row_end && Aj[jj] == j ){ + x += Ax[jj]; + jj++; + } + Aj[nnz] = j; + Ax[nnz] = x; + nnz++; + } + Ap[i+1] = nnz; + } +} + +EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals) +{ + using namespace Eigen; + SparseMatrix<Scalar> mat(SIZE,SIZE); + mat.resizeNonZeros(coords.size()); +// std::cerr << "setrand_scipy...\n"; + coo_tocsr<Scalar>(SIZE,SIZE, coords.size(), coords, vals, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr()); +// std::cerr << "coo_tocsr ok\n"; + + csr_sort_indices(SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr()); + + csr_sum_duplicates(SIZE, SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr()); + + mat.resizeNonZeros(mat._outerIndexPtr()[SIZE]); + + return &mat.coeffRef(coords[0].x(), coords[0].y()); +} + + #ifndef NOUBLAS EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals) { @@ -242,7 +430,7 @@ mapped_matrix<Scalar> aux(SIZE,SIZE); for (int i=0; i<coords.size(); ++i) { - aux(coords[i].x(), coords[i].y()) = vals[i]; + aux(coords[i].x(), coords[i].y()) += vals[i]; } CHECK_MEM; compressed_matrix<Scalar> mat(aux); @@ -278,12 +466,12 @@ using namespace boost; using namespace boost::numeric; using namespace boost::numeric::ublas; - + // ublas::vector<coordinate_vector<Scalar> > foo; generalized_vector_of_vector<Scalar, row_major, ublas::vector<coordinate_vector<Scalar> > > aux(SIZE,SIZE); for (int i=0; i<coords.size(); ++i) { - aux(coords[i].x(), coords[i].y()) = vals[i]; + aux(coords[i].x(), coords[i].y()) += vals[i]; } CHECK_MEM; compressed_matrix<Scalar,row_major> mat(aux);
diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt new file mode 100644 index 0000000..a6c330a --- /dev/null +++ b/blas/CMakeLists.txt
@@ -0,0 +1,14 @@ +project(EigenBlas) + +add_custom_target(blas) + +set(EigenBlas_SRCS single.cpp double.cpp complex_single.cpp complex_double.cpp) + +add_library(eigen_blas SHARED ${EigenBlas_SRCS}) +add_dependencies(blas eigen_blas) + +install(TARGETS eigen_blas + RUNTIME DESTINATION bin + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib) +
diff --git a/blas/README.txt b/blas/README.txt new file mode 100644 index 0000000..466a675 --- /dev/null +++ b/blas/README.txt
@@ -0,0 +1,7 @@ + +This directory contains a BLAS library built on top of Eigen. + +This is currently a work in progress which is far to be ready for use, +but feel free to contribute to it if you wish. + +If you want to compile it, set the cmake variable EIGEN_BUILD_BLAS to "on".
diff --git a/blas/common.h b/blas/common.h new file mode 100644 index 0000000..74c3c9f --- /dev/null +++ b/blas/common.h
@@ -0,0 +1,115 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_BLAS_COMMON_H +#define EIGEN_BLAS_COMMON_H + +#ifndef SCALAR +#error the token SCALAR must be defined to compile this file +#endif + +#ifdef __cplusplus +extern "C" +{ +#endif + +#include <blas.h> + +#ifdef __cplusplus +} +#endif + + +#define NOTR 0 +#define TR 1 +#define ADJ 2 + +#define LEFT 0 +#define RIGHT 1 + +#define UP 0 +#define LO 1 + +#define NUNIT 0 +#define UNIT 1 + +#define OP(X) ( ((X)=='N' || (X)=='n') ? NOTR \ + : ((X)=='T' || (X)=='t') ? TR \ + : ((X)=='C' || (X)=='c') ? ADJ \ + : 0xff) + +#define SIDE(X) ( ((X)=='L' || (X)=='l') ? LEFT \ + : ((X)=='R' || (X)=='r') ? RIGHT \ + : 0xff) + +#define UPLO(X) ( ((X)=='U' || (X)=='u') ? UP \ + : ((X)=='L' || (X)=='l') ? LO \ + : 0xff) + +#define DIAG(X) ( ((X)=='N' || (X)=='N') ? NUNIT \ + : ((X)=='U' || (X)=='u') ? UNIT \ + : 0xff) + +#include <Eigen/Core> +#include <Eigen/Jacobi> +using namespace Eigen; + +template<typename T> +Block<NestByValue<Map<Matrix<T,Dynamic,Dynamic> > >, Dynamic, Dynamic> +matrix(T* data, int rows, int cols, int stride) +{ + return Map<Matrix<T,Dynamic,Dynamic> >(data, stride, cols).nestByValue().block(0,0,rows,cols); +} + +template<typename T> +Block<NestByValue<Map<Matrix<T,Dynamic,Dynamic,RowMajor> > >, Dynamic, 1> +vector(T* data, int size, int incr) +{ + return Map<Matrix<T,Dynamic,Dynamic,RowMajor> >(data, size, incr).nestByValue().col(0); +} + +template<typename T> +Map<Matrix<T,Dynamic,1> > +vector(T* data, int size) +{ + return Map<Matrix<T,Dynamic,1> >(data, size); +} + +typedef SCALAR Scalar; +typedef NumTraits<Scalar>::Real RealScalar; +typedef std::complex<RealScalar> Complex; + +enum +{ + IsComplex = Eigen::NumTraits<SCALAR>::IsComplex, + Conj = IsComplex +}; + +typedef Block<NestByValue<Map<Matrix<Scalar,Dynamic,Dynamic> > >, Dynamic, Dynamic> MatrixType; +typedef Block<NestByValue<Map<Matrix<Scalar,Dynamic,Dynamic, RowMajor> > >, Dynamic, 1> StridedVectorType; +typedef Map<Matrix<Scalar,Dynamic,1> > CompactVectorType; + +#define EIGEN_BLAS_FUNC(X) EIGEN_CAT(SCALAR_SUFFIX,X##_) + +#endif // EIGEN_BLAS_COMMON_H
diff --git a/Eigen/src/QR/QrInstantiations.cpp b/blas/complex_double.cpp similarity index 79% copy from Eigen/src/QR/QrInstantiations.cpp copy to blas/complex_double.cpp index 695377d..f51ccb2 100644 --- a/Eigen/src/QR/QrInstantiations.cpp +++ b/blas/complex_double.cpp
@@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -22,17 +22,10 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifndef EIGEN_EXTERN_INSTANTIATIONS -#define EIGEN_EXTERN_INSTANTIATIONS -#endif -#include "../../Core" -#undef EIGEN_EXTERN_INSTANTIATIONS +#define SCALAR std::complex<double> +#define SCALAR_SUFFIX c +#define ISCOMPLEX 1 -#include "../../QR" - -namespace Eigen -{ - -EIGEN_QR_MODULE_INSTANTIATE(); - -} +#include "level1_impl.h" +#include "level2_impl.h" +#include "level3_impl.h"
diff --git a/Eigen/src/QR/QrInstantiations.cpp b/blas/complex_single.cpp similarity index 79% copy from Eigen/src/QR/QrInstantiations.cpp copy to blas/complex_single.cpp index 695377d..b6617e7 100644 --- a/Eigen/src/QR/QrInstantiations.cpp +++ b/blas/complex_single.cpp
@@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -22,17 +22,10 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifndef EIGEN_EXTERN_INSTANTIATIONS -#define EIGEN_EXTERN_INSTANTIATIONS -#endif -#include "../../Core" -#undef EIGEN_EXTERN_INSTANTIATIONS +#define SCALAR std::complex<float> +#define SCALAR_SUFFIX z +#define ISCOMPLEX 1 -#include "../../QR" - -namespace Eigen -{ - -EIGEN_QR_MODULE_INSTANTIATE(); - -} +#include "level1_impl.h" +#include "level2_impl.h" +#include "level3_impl.h"
diff --git a/Eigen/src/QR/QrInstantiations.cpp b/blas/double.cpp similarity index 79% copy from Eigen/src/QR/QrInstantiations.cpp copy to blas/double.cpp index 695377d..8145696 100644 --- a/Eigen/src/QR/QrInstantiations.cpp +++ b/blas/double.cpp
@@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -22,17 +22,10 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifndef EIGEN_EXTERN_INSTANTIATIONS -#define EIGEN_EXTERN_INSTANTIATIONS -#endif -#include "../../Core" -#undef EIGEN_EXTERN_INSTANTIATIONS +#define SCALAR double +#define SCALAR_SUFFIX d +#define ISCOMPLEX 0 -#include "../../QR" - -namespace Eigen -{ - -EIGEN_QR_MODULE_INSTANTIATE(); - -} +#include "level1_impl.h" +#include "level2_impl.h" +#include "level3_impl.h"
diff --git a/blas/level1_impl.h b/blas/level1_impl.h new file mode 100644 index 0000000..c508626 --- /dev/null +++ b/blas/level1_impl.h
@@ -0,0 +1,225 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#include "common.h" + +int EIGEN_BLAS_FUNC(axpy)(int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + if(*incx==1 && *incy==1) + vector(y,*n) += alpha * vector(x,*n); + else + vector(y,*n,*incy) += alpha * vector(x,*n,*incx); + + return 1; +} + +// computes the sum of magnitudes of all vector elements or, for a complex vector x, the sum +// res = |Rex1| + |Imx1| + |Rex2| + |Imx2| + ... + |Rexn| + |Imxn|, where x is a vector of order n +RealScalar EIGEN_BLAS_FUNC(asum)(int *n, RealScalar *px, int *incx) +{ + int size = IsComplex ? 2* *n : *n; + + if(*incx==1) + return vector(px,size).cwise().abs().sum(); + else + return vector(px,size,*incx).cwise().abs().sum(); + + return 1; +} + +int EIGEN_BLAS_FUNC(copy)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ + int size = IsComplex ? 2* *n : *n; + + if(*incx==1 && *incy==1) + vector(py,size) = vector(px,size); + else + vector(py,size,*incy) = vector(px,size,*incx); + + return 1; +} + +// computes a vector-vector dot product. +Scalar EIGEN_BLAS_FUNC(dot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + + if(*incx==1 && *incy==1) + return (vector(x,*n).cwise()*vector(y,*n)).sum(); + + return (vector(x,*n,*incx).cwise()*vector(y,*n,*incy)).sum(); +} + +/* + +// computes a vector-vector dot product with extended precision. +Scalar EIGEN_BLAS_FUNC(sdot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ + // TODO + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + + if(*incx==1 && *incy==1) + return vector(x,*n).dot(vector(y,*n)); + + return vector(x,*n,*incx).dot(vector(y,*n,*incy)); +} + +*/ + +#if ISCOMPLEX + +// computes a dot product of a conjugated vector with another vector. +Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + + if(*incx==1 && *incy==1) + return vector(x,*n).dot(vector(y,*n)); + + return vector(x,*n,*incx).dot(vector(y,*n,*incy)); +} + +// computes a vector-vector dot product without complex conjugation. +Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + + if(*incx==1 && *incy==1) + return (vector(x,*n).cwise()*vector(y,*n)).sum(); + + return (vector(x,*n,*incx).cwise()*vector(y,*n,*incy)).sum(); +} + +#endif // ISCOMPLEX + +// computes the Euclidean norm of a vector. +Scalar EIGEN_BLAS_FUNC(nrm2)(int *n, RealScalar *px, int *incx) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + + if(*incx==1) + return vector(x,*n).norm(); + + return vector(x,*n,*incx).norm(); +} + +int EIGEN_BLAS_FUNC(rot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + Scalar c = *reinterpret_cast<Scalar*>(pc); + Scalar s = *reinterpret_cast<Scalar*>(ps); + + StridedVectorType vx(vector(x,*n,*incx)); + StridedVectorType vy(vector(y,*n,*incy)); + ei_apply_rotation_in_the_plane(vx, vy, PlanarRotation<Scalar>(c,s)); + return 1; +} + +int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps) +{ + Scalar a = *reinterpret_cast<Scalar*>(pa); + Scalar b = *reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar* s = reinterpret_cast<Scalar*>(ps); + + PlanarRotation<Scalar> r; + r.makeGivens(a,b); + *c = r.c(); + *s = r.s(); + + return 1; +} + +#if !ISCOMPLEX +/* +// performs rotation of points in the modified plane. +int EIGEN_BLAS_FUNC(rotm)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *param) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + + // TODO + + return 0; +} + +// computes the modified parameters for a Givens rotation. +int EIGEN_BLAS_FUNC(rotmg)(RealScalar *d1, RealScalar *d2, RealScalar *x1, RealScalar *x2, RealScalar *param) +{ + // TODO + + return 0; +} +*/ +#endif // !ISCOMPLEX + +int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *px, int *incx, RealScalar *palpha) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + if(*incx==1) + vector(x,*n) *= alpha; + + vector(x,*n,*incx) *= alpha; + + return 1; +} + +int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ + int size = IsComplex ? 2* *n : *n; + + if(*incx==1 && *incy==1) + vector(py,size).swap(vector(px,size)); + else + vector(py,size,*incy).swap(vector(px,size,*incx)); + + return 1; +} + +#if !ISCOMPLEX + +RealScalar EIGEN_BLAS_FUNC(casum)(int *n, RealScalar *px, int *incx) +{ + Complex* x = reinterpret_cast<Complex*>(px); + + if(*incx==1) + return vector(x,*n).cwise().abs().sum(); + else + return vector(x,*n,*incx).cwise().abs().sum(); + + return 1; +} + +#endif // ISCOMPLEX
diff --git a/blas/level2_impl.h b/blas/level2_impl.h new file mode 100644 index 0000000..5691e8a --- /dev/null +++ b/blas/level2_impl.h
@@ -0,0 +1,214 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#include "common.h" + +int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc) +{ + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + if(beta!=Scalar(1)) + vector(c, *m, *incc) *= beta; + + if(OP(*opa)==NOTR) + if(*incc==1) + vector(c,*m) += alpha * matrix(a,*m,*n,*lda) * vector(b,*n,*incb); + else + vector(c,*m,*incc) += alpha * matrix(a,*m,*n,*lda) * vector(b,*n,*incb); + else if(OP(*opa)==TR) + if(*incb==1) + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).transpose() * vector(b,*n); + else + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).transpose() * vector(b,*n,*incb); + else if(OP(*opa)==TR) + if(*incb==1) + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).adjoint() * vector(b,*n); + else + vector(c,*m,*incc) += alpha * matrix(a,*n,*m,*lda).adjoint() * vector(b,*n,*incb); + else + return 0; + + return 1; +} + +/* +int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb) +{ + typedef void (*functype)(int, const Scalar *, int, Scalar *, int); + functype func[16]; + + static bool init = false; + if(!init) + { + for(int k=0; k<16; ++k) + func[k] = 0; + +// func[NOTR | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|0, false,ColMajor,ColMajor>::run); +// func[TR | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|0, false,RowMajor,ColMajor>::run); +// func[ADJ | (UP << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|0, Conj, RowMajor,ColMajor>::run); +// +// func[NOTR | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|0, false,ColMajor,ColMajor>::run); +// func[TR | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|0, false,RowMajor,ColMajor>::run); +// func[ADJ | (LO << 2) | (NUNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|0, Conj, RowMajor,ColMajor>::run); +// +// func[NOTR | (UP << 3) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run); +// func[TR | (UP << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run); +// func[ADJ | (UP << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, UpperTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run); +// +// func[NOTR | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run); +// func[TR | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run); +// func[ADJ | (LO << 2) | (UNIT << 3)] = (ei_triangular_solve_vector<Scalar, LowerTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + + int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); + if(code>=16 || func[code]==0) + return 0; + + func[code](*n, a, *lda, b, *incb); + return 1; +} +*/ + +/* +int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb) +{ + // TODO + + typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int); + functype func[16]; + + static bool init = false; + if(!init) + { + for(int k=0; k<16; ++k) + func[k] = 0; + +// func[NOTR | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run); +// func[TR | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run); +// func[ADJ | (UP << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run); +// +// func[NOTR | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run); +// func[TR | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run); +// func[ADJ | (LO << 2) | (NUNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run); +// +// func[NOTR | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run); +// func[TR | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run); +// func[ADJ | (UP << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run); +// +// func[NOTR | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run); +// func[TR | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run); +// func[ADJ | (LO << 2) | (UNIT << 3)] = (ei_product_triangular_matrix_vector<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + + int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); + if(code>=16 || func[code]==0) + return 0; + + func[code](*n, a, *lda, b, *incb, b, *incb); + return 1; +} +*/ + +/* +int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *inca, RealScalar *pc, int *ldc) +{ + // TODO + typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar); + functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + +// func[UP] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run); +// func[LO] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, a, *inca, c, *ldc, alpha); + return 1; +} +*/ + +/* +int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *inca, RealScalar *pb, int *incb, RealScalar *pc, int *ldc) +{ + // TODO + typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar); + functype func[2]; + + static bool init = false; + if(!init) + { + for(int k=0; k<2; ++k) + func[k] = 0; + +// func[UP] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run); +// func[LO] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + int code = UPLO(*uplo); + if(code>=2 || func[code]==0) + return 0; + + func[code](*n, a, *inca, b, *incb, c, *ldc, alpha); + return 1; +} +*/ + +#if ISCOMPLEX + +#endif // ISCOMPLEX
diff --git a/blas/level3_impl.h b/blas/level3_impl.h new file mode 100644 index 0000000..d44de1b --- /dev/null +++ b/blas/level3_impl.h
@@ -0,0 +1,365 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#include "common.h" + +int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + typedef void (*functype)(int, int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar); + functype func[12]; + + static bool init = false; + if(!init) + { + for(int k=0; k<12; ++k) + func[k] = 0; + func[NOTR | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,ColMajor,false,ColMajor>::run); + func[ADJ | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, ColMajor,false,ColMajor>::run); + func[NOTR | (TR << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,RowMajor,false,ColMajor>::run); + func[TR | (TR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,RowMajor,false,ColMajor>::run); + func[ADJ | (TR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, RowMajor,false,ColMajor>::run); + func[NOTR | (ADJ << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,RowMajor,Conj, ColMajor>::run); + func[TR | (ADJ << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,RowMajor,Conj, ColMajor>::run); + func[ADJ | (ADJ << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, RowMajor,Conj, ColMajor>::run); + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + if(beta!=Scalar(1)) + matrix(c, *m, *n, *ldc) *= beta; + + int code = OP(*opa) | (OP(*opb) << 2); + if(code>=12 || func[code]==0) + return 0; + + func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha); + return 1; +} + +int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb) +{ + typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int); + functype func[32]; + + static bool init = false; + if(!init) + { + for(int k=0; k<32; ++k) + func[k] = 0; + + func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|0, false,ColMajor,ColMajor>::run); + func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|0, false,RowMajor,ColMajor>::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|0, Conj, RowMajor,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|0, false,ColMajor,ColMajor>::run); + func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|0, false,RowMajor,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|0, Conj, RowMajor,ColMajor>::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|0, false,ColMajor,ColMajor>::run); + func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|0, false,RowMajor,ColMajor>::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|0, Conj, RowMajor,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|0, false,ColMajor,ColMajor>::run); + func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|0, false,RowMajor,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|0, Conj, RowMajor,ColMajor>::run); + + + func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run); + func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, UpperTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run); + func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,UpperTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run); + func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, LowerTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|UnitDiagBit,false,ColMajor,ColMajor>::run); + func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|UnitDiagBit,false,RowMajor,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,LowerTriangular|UnitDiagBit,Conj, RowMajor,ColMajor>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + // TODO handle alpha + + int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4); + if(code>=32 || func[code]==0) + return 0; + + func[code](*m, *n, a, *lda, b, *ldb); + return 1; +} + + +// b = alpha*op(a)*b for side = 'L'or'l' +// b = alpha*b*op(a) for side = 'R'or'r' +int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb) +{ + typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar); + functype func[32]; + + static bool init = false; + if(!init) + { + for(int k=0; k<32; ++k) + func[k] = 0; + + func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, false,ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, false,ColMajor,false,RowMajor,false,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, true, ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, true, RowMajor,false,ColMajor,false,ColMajor>::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, false,ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, false,ColMajor,false,RowMajor,false,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run); + + func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run); + func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,false,ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,false,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,UpperTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,Conj, ColMajor>::run); + + func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,true, ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,false,ColMajor,false,ColMajor>::run); + func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,true, RowMajor,Conj, ColMajor,false,ColMajor>::run); + + func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,false,ColMajor,false,ColMajor,false,ColMajor>::run); + func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,false,ColMajor>::run); + func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,LowerTriangular|UnitDiagBit,false,ColMajor,false,RowMajor,Conj, ColMajor>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + + int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4); + if(code>=32 || func[code]==0) + return 0; + + func[code](*m, *n, a, *lda, b, *ldb, b, *ldb, alpha); + return 1; +} + +// c = alpha*a*b + beta*c for side = 'L'or'l' +// c = alpha*b*a + beta*c for side = 'R'or'r +int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + if(beta!=Scalar(1)) + matrix(c, *m, *n, *ldc) *= beta; + + if(SIDE(*side)==LEFT) + if(UPLO(*uplo)==UP) + ei_product_selfadjoint_matrix<Scalar, RowMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else if(UPLO(*uplo)==LO) + ei_product_selfadjoint_matrix<Scalar, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else + return 0; + else if(SIDE(*side)==RIGHT) + if(UPLO(*uplo)==UP) + ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, RowMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); + else if(UPLO(*uplo)==LO) + ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); + else + return 0; + else + return 0; + + return 1; +} + +// c = alpha*a*a' + beta*c for op = 'N'or'n' +// c = alpha*a'*a + beta*c for op = 'T'or't','C'or'c' +int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int, Scalar); + functype func[8]; + + static bool init = false; + if(!init) + { + for(int k=0; k<8; ++k) + func[k] = 0; + + func[NOTR | (UP << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, UpperTriangular>::run); + func[TR | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,UpperTriangular>::run); + func[ADJ | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,UpperTriangular>::run); + + func[NOTR | (LO << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, LowerTriangular>::run); + func[TR | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,LowerTriangular>::run); + func[ADJ | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,LowerTriangular>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + int code = OP(*op) | (UPLO(*uplo) << 2); + if(code>=8 || func[code]==0) + return 0; + + if(beta!=Scalar(1)) + matrix(c, *n, *n, *ldc) *= beta; + + func[code](*n, *k, a, *lda, c, *ldc, alpha); + return 1; +} + +// c = alpha*a*b' + alpha*b*a' + beta*c for op = 'N'or'n' +// c = alpha*a'*b + alpha*b'*a + beta*c for op = 'T'or't' +int EIGEN_BLAS_FUNC(syr2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + // TODO + + return 0; +} + + +#if ISCOMPLEX + +// c = alpha*a*b + beta*c for side = 'L'or'l' +// c = alpha*b*a + beta*c for side = 'R'or'r +int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + if(beta!=Scalar(1)) + matrix(c, *m, *n, *ldc) *= beta; + + if(SIDE(*side)==LEFT) + if(UPLO(*uplo)==UP) + ei_product_selfadjoint_matrix<Scalar, RowMajor,true,Conj, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else if(UPLO(*uplo)==LO) + ei_product_selfadjoint_matrix<Scalar, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); + else + return 0; + else if(SIDE(*side)==RIGHT) + if(UPLO(*uplo)==UP) + ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, RowMajor,true,Conj, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); + else if(UPLO(*uplo)==LO) + ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); + else + return 0; + else + return 0; + + return 1; +} + +// c = alpha*a*conj(a') + beta*c for op = 'N'or'n' +// c = alpha*conj(a')*a + beta*c for op = 'C'or'c' +int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int, Scalar); + functype func[8]; + + static bool init = false; + if(!init) + { + for(int k=0; k<8; ++k) + func[k] = 0; + + func[NOTR | (UP << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, UpperTriangular>::run); + func[ADJ | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,UpperTriangular>::run); + + func[NOTR | (LO << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, LowerTriangular>::run); + func[ADJ | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,LowerTriangular>::run); + + init = true; + } + + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + int code = OP(*op) | (UPLO(*uplo) << 2); + if(code>=8 || func[code]==0) + return 0; + + if(beta!=Scalar(1)) + matrix(c, *n, *n, *ldc) *= beta; + + func[code](*n, *k, a, *lda, c, *ldc, alpha); + return 1; +} + +// c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c, for op = 'N'or'n' +// c = alpha*conj(b')*a + conj(alpha)*conj(a')*b + beta*c, for op = 'C'or'c' +int EIGEN_BLAS_FUNC(her2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) +{ + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* b = reinterpret_cast<Scalar*>(pb); + Scalar* c = reinterpret_cast<Scalar*>(pc); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + // TODO + + return 0; +} + +#endif // ISCOMPLEX
diff --git a/Eigen/src/QR/QrInstantiations.cpp b/blas/single.cpp similarity index 79% rename from Eigen/src/QR/QrInstantiations.cpp rename to blas/single.cpp index 695377d..842e104 100644 --- a/Eigen/src/QR/QrInstantiations.cpp +++ b/blas/single.cpp
@@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -22,17 +22,10 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifndef EIGEN_EXTERN_INSTANTIATIONS -#define EIGEN_EXTERN_INSTANTIATIONS -#endif -#include "../../Core" -#undef EIGEN_EXTERN_INSTANTIATIONS +#define SCALAR float +#define SCALAR_SUFFIX s +#define ISCOMPLEX 0 -#include "../../QR" - -namespace Eigen -{ - -EIGEN_QR_MODULE_INSTANTIATE(); - -} +#include "level1_impl.h" +#include "level2_impl.h" +#include "level3_impl.h"
diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake index faa75c6..b8e159b 100644 --- a/cmake/EigenTesting.cmake +++ b/cmake/EigenTesting.cmake
@@ -27,7 +27,11 @@ # void test_<testname>() { ... } # # this macro add an executable test_<testname> as well as a ctest test -# named <testname> +# named <testname>. +# +# it also adds another executable debug_<testname> that compiles in full debug mode +# and is not added to the test target. The idea is that when a test fails you want +# a quick way of rebuilding this specific test in full debug mode. # # On platforms with bash simply run: # "ctest -V" or "ctest -V -R <testname>" @@ -36,38 +40,65 @@ macro(ei_add_test testname) set(targetname test_${testname}) + if(NOT MSVC_IDE) + set(debug_targetname debug_${testname}) + endif(NOT MSVC_IDE) set(filename ${testname}.cpp) add_executable(${targetname} ${filename}) + add_dependencies(btest ${targetname}) + if(NOT MSVC_IDE) + add_executable(${debug_targetname} EXCLUDE_FROM_ALL ${filename}) + endif(NOT MSVC_IDE) if(NOT EIGEN_NO_ASSERTION_CHECKING) if(MSVC) set_target_properties(${targetname} PROPERTIES COMPILE_FLAGS "/EHsc") + if(NOT MSVC_IDE) + set_target_properties(${debug_targetname} PROPERTIES COMPILE_FLAGS "/EHsc") + endif(NOT MSVC_IDE) else(MSVC) set_target_properties(${targetname} PROPERTIES COMPILE_FLAGS "-fexceptions") + set_target_properties(${debug_targetname} PROPERTIES COMPILE_FLAGS "-fexceptions") endif(MSVC) option(EIGEN_DEBUG_ASSERTS "Enable debuging of assertions" OFF) if(EIGEN_DEBUG_ASSERTS) ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_DEBUG_ASSERTS=1") + if(NOT MSVC_IDE) + ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-DEIGEN_DEBUG_ASSERTS=1") + endif(NOT MSVC_IDE) endif(EIGEN_DEBUG_ASSERTS) else(NOT EIGEN_NO_ASSERTION_CHECKING) ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_NO_ASSERTION_CHECKING=1") + if(NOT MSVC_IDE) + ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-DEIGEN_NO_ASSERTION_CHECKING=1") + endif(NOT MSVC_IDE) endif(NOT EIGEN_NO_ASSERTION_CHECKING) + # let the user pass e.g. optimization flags, but don't apply them to the debug target if(${ARGC} GREATER 1) ei_add_target_property(${targetname} COMPILE_FLAGS "${ARGV1}") endif(${ARGC} GREATER 1) - ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}") + # for the debug target, add full debug options + if(CMAKE_COMPILER_IS_GNUCXX) + # O0 is in principle redundant here, but doesn't hurt + ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-O0 -g3") + elseif(MSVC) + if(NOT MSVC_IDE) + ei_add_target_property(${debug_targetname} COMPILE_FLAGS "/Od /Zi") + endif(NOT MSVC_IDE) + endif(CMAKE_COMPILER_IS_GNUCXX) - if(TEST_LIB) - target_link_libraries(${targetname} Eigen2) - endif(TEST_LIB) + ei_add_target_property(${targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}") + if(NOT MSVC_IDE) + ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-DEIGEN_TEST_FUNC=${testname}") + endif(NOT MSVC_IDE) target_link_libraries(${targetname} ${EXTERNAL_LIBS}) if(${ARGC} GREATER 2) @@ -75,11 +106,18 @@ string(LENGTH "${ARGV2_stripped}" ARGV2_stripped_length) if(${ARGV2_stripped_length} GREATER 0) target_link_libraries(${targetname} ${ARGV2}) + if(NOT MSVC_IDE) + target_link_libraries(${debug_targetname} ${ARGV2}) + endif(NOT MSVC_IDE) endif(${ARGV2_stripped_length} GREATER 0) endif(${ARGC} GREATER 2) if(WIN32) - add_test(${testname} "${targetname}") + if(CYGWIN) + add_test(${testname} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname}") + else(CYGWIN) + add_test(${testname} "${targetname}") + endif(CYGWIN) else(WIN32) add_test(${testname} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname}") endif(WIN32) @@ -102,31 +140,31 @@ if(EIGEN_TEST_SSE2) message("SSE2: ON") else(EIGEN_TEST_SSE2) - message("SSE2: AUTO") + message("SSE2: Using architecture defaults") endif(EIGEN_TEST_SSE2) if(EIGEN_TEST_SSE3) message("SSE3: ON") else(EIGEN_TEST_SSE3) - message("SSE3: AUTO") + message("SSE3: Using architecture defaults") endif(EIGEN_TEST_SSE3) if(EIGEN_TEST_SSSE3) message("SSSE3: ON") else(EIGEN_TEST_SSSE3) - message("SSSE3: AUTO") + message("SSSE3: Using architecture defaults") endif(EIGEN_TEST_SSSE3) if(EIGEN_TEST_ALTIVEC) - message("Altivec: ON") + message("Altivec: Using architecture defaults") else(EIGEN_TEST_ALTIVEC) - message("Altivec: AUTO") + message("Altivec: Using architecture defaults") endif(EIGEN_TEST_ALTIVEC) if(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION) message("Explicit vec: OFF") else(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION) - message("Explicit vec: AUTO") + message("Explicit vec: Using architecture defaults") endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION) message("\n${EIGEN_TESTING_SUMMARY}") @@ -159,6 +197,9 @@ else(EIGEN_COVERAGE_TESTING) set(COVERAGE_FLAGS "") endif(EIGEN_COVERAGE_TESTING) + if(EIGEN_TEST_RVALUE_REF_SUPPORT OR EIGEN_TEST_C++0x) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x") + endif(EIGEN_TEST_RVALUE_REF_SUPPORT OR EIGEN_TEST_C++0x) if(CMAKE_SYSTEM_NAME MATCHES Linux) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -g2") set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} ${COVERAGE_FLAGS} -O2 -g2")
diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index 4e8f516..ce25125 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt
@@ -1,3 +1,11 @@ +project(EigenDemos) -add_subdirectory(mandelbrot) -add_subdirectory(opengl) +add_custom_target(demos) + +find_package(Qt4) +if(QT4_FOUND) + add_subdirectory(mandelbrot) + add_subdirectory(opengl) +else(QT4_FOUND) + message(STATUS "Qt4 not found, so disabling the mandelbrot and opengl demos") +endif(QT4_FOUND)
diff --git a/demos/mandelbrot/CMakeLists.txt b/demos/mandelbrot/CMakeLists.txt index d34b60a..5c500e0 100644 --- a/demos/mandelbrot/CMakeLists.txt +++ b/demos/mandelbrot/CMakeLists.txt
@@ -16,5 +16,6 @@ qt4_automoc(${mandelbrot_SRCS}) add_executable(mandelbrot ${mandelbrot_SRCS}) +add_dependencies(demos mandelbrot) target_link_libraries(mandelbrot ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY})
diff --git a/demos/mix_eigen_and_c/README b/demos/mix_eigen_and_c/README new file mode 100644 index 0000000..21dba86 --- /dev/null +++ b/demos/mix_eigen_and_c/README
@@ -0,0 +1,9 @@ +This is an example of how one can wrap some of Eigen into a C library. + +To try this with GCC, do: + + g++ -c binary_library.cpp -O2 -msse2 -I ../.. + gcc example.c binary_library.o -o example -lstdc++ + ./example + +TODO: add CMakeLists, add more explanations here \ No newline at end of file
diff --git a/demos/mix_eigen_and_c/binary_library.cpp b/demos/mix_eigen_and_c/binary_library.cpp new file mode 100644 index 0000000..efbf28f --- /dev/null +++ b/demos/mix_eigen_and_c/binary_library.cpp
@@ -0,0 +1,200 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <http://www.gnu.org/licenses/>. + +// This C++ file compiles to binary code that can be linked to by your C program, +// thanks to the extern "C" syntax used in the declarations in binary_library.h. + +#include "binary_library.h" + +#include <Eigen/Core> + +using namespace Eigen; + +/************************* pointer conversion methods **********************************************/ + +////// class MatrixXd ////// + +inline MatrixXd& c_to_eigen(C_MatrixXd* ptr) +{ + return *reinterpret_cast<MatrixXd*>(ptr); +} + +inline const MatrixXd& c_to_eigen(const C_MatrixXd* ptr) +{ + return *reinterpret_cast<const MatrixXd*>(ptr); +} + +inline C_MatrixXd* eigen_to_c(MatrixXd& ref) +{ + return reinterpret_cast<C_MatrixXd*>(&ref); +} + +inline const C_MatrixXd* eigen_to_c(const MatrixXd& ref) +{ + return reinterpret_cast<const C_MatrixXd*>(&ref); +} + +////// class Map<MatrixXd> ////// + +inline Map<MatrixXd>& c_to_eigen(C_Map_MatrixXd* ptr) +{ + return *reinterpret_cast<Map<MatrixXd>*>(ptr); +} + +inline const Map<MatrixXd>& c_to_eigen(const C_Map_MatrixXd* ptr) +{ + return *reinterpret_cast<const Map<MatrixXd>*>(ptr); +} + +inline C_Map_MatrixXd* eigen_to_c(Map<MatrixXd>& ref) +{ + return reinterpret_cast<C_Map_MatrixXd*>(&ref); +} + +inline const C_Map_MatrixXd* eigen_to_c(const Map<MatrixXd>& ref) +{ + return reinterpret_cast<const C_Map_MatrixXd*>(&ref); +} + + +/************************* implementation of classes **********************************************/ + + +////// class MatrixXd ////// + + +C_MatrixXd* MatrixXd_new(int rows, int cols) +{ + return eigen_to_c(*new MatrixXd(rows,cols)); +} + +void MatrixXd_delete(C_MatrixXd *m) +{ + delete &c_to_eigen(m); +} + +double* MatrixXd_data(C_MatrixXd *m) +{ + return c_to_eigen(m).data(); +} + +void MatrixXd_set_zero(C_MatrixXd *m) +{ + c_to_eigen(m).setZero(); +} + +void MatrixXd_resize(C_MatrixXd *m, int rows, int cols) +{ + c_to_eigen(m).resize(rows,cols); +} + +void MatrixXd_copy(C_MatrixXd *dst, const C_MatrixXd *src) +{ + c_to_eigen(dst) = c_to_eigen(src); +} + +void MatrixXd_copy_map(C_MatrixXd *dst, const C_Map_MatrixXd *src) +{ + c_to_eigen(dst) = c_to_eigen(src); +} + +void MatrixXd_set_coeff(C_MatrixXd *m, int i, int j, double coeff) +{ + c_to_eigen(m)(i,j) = coeff; +} + +double MatrixXd_get_coeff(const C_MatrixXd *m, int i, int j) +{ + return c_to_eigen(m)(i,j); +} + +void MatrixXd_print(const C_MatrixXd *m) +{ + std::cout << c_to_eigen(m) << std::endl; +} + +void MatrixXd_multiply(const C_MatrixXd *m1, const C_MatrixXd *m2, C_MatrixXd *result) +{ + c_to_eigen(result) = c_to_eigen(m1) * c_to_eigen(m2); +} + +void MatrixXd_add(const C_MatrixXd *m1, const C_MatrixXd *m2, C_MatrixXd *result) +{ + c_to_eigen(result) = c_to_eigen(m1) + c_to_eigen(m2); +} + + + +////// class Map_MatrixXd ////// + + +C_Map_MatrixXd* Map_MatrixXd_new(double *array, int rows, int cols) +{ + return eigen_to_c(*new Map<MatrixXd>(array,rows,cols)); +} + +void Map_MatrixXd_delete(C_Map_MatrixXd *m) +{ + delete &c_to_eigen(m); +} + +void Map_MatrixXd_set_zero(C_Map_MatrixXd *m) +{ + c_to_eigen(m).setZero(); +} + +void Map_MatrixXd_copy(C_Map_MatrixXd *dst, const C_Map_MatrixXd *src) +{ + c_to_eigen(dst) = c_to_eigen(src); +} + +void Map_MatrixXd_copy_matrix(C_Map_MatrixXd *dst, const C_MatrixXd *src) +{ + c_to_eigen(dst) = c_to_eigen(src); +} + +void Map_MatrixXd_set_coeff(C_Map_MatrixXd *m, int i, int j, double coeff) +{ + c_to_eigen(m)(i,j) = coeff; +} + +double Map_MatrixXd_get_coeff(const C_Map_MatrixXd *m, int i, int j) +{ + return c_to_eigen(m)(i,j); +} + +void Map_MatrixXd_print(const C_Map_MatrixXd *m) +{ + std::cout << c_to_eigen(m) << std::endl; +} + +void Map_MatrixXd_multiply(const C_Map_MatrixXd *m1, const C_Map_MatrixXd *m2, C_Map_MatrixXd *result) +{ + c_to_eigen(result) = c_to_eigen(m1) * c_to_eigen(m2); +} + +void Map_MatrixXd_add(const C_Map_MatrixXd *m1, const C_Map_MatrixXd *m2, C_Map_MatrixXd *result) +{ + c_to_eigen(result) = c_to_eigen(m1) + c_to_eigen(m2); +}
diff --git a/demos/mix_eigen_and_c/binary_library.h b/demos/mix_eigen_and_c/binary_library.h new file mode 100644 index 0000000..3a4dc5b --- /dev/null +++ b/demos/mix_eigen_and_c/binary_library.h
@@ -0,0 +1,86 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <http://www.gnu.org/licenses/>. + +// This is a pure C header, no C++ here. +// The functions declared here will be implemented in C++ but +// we don't have to know, because thanks to the extern "C" syntax, +// they will be compiled to C object code. + +#ifdef __cplusplus +extern "C" +{ +#endif + + // just dummy empty structs to give different pointer types, + // instead of using void* which would be type unsafe + struct C_MatrixXd {}; + struct C_Map_MatrixXd {}; + + // the C_MatrixXd class, wraps some of the functionality + // of Eigen::MatrixXd. + struct C_MatrixXd* MatrixXd_new(int rows, int cols); + void MatrixXd_delete (struct C_MatrixXd *m); + double* MatrixXd_data (struct C_MatrixXd *m); + void MatrixXd_set_zero (struct C_MatrixXd *m); + void MatrixXd_resize (struct C_MatrixXd *m, int rows, int cols); + void MatrixXd_copy (struct C_MatrixXd *dst, + const struct C_MatrixXd *src); + void MatrixXd_copy_map (struct C_MatrixXd *dst, + const struct C_Map_MatrixXd *src); + void MatrixXd_set_coeff (struct C_MatrixXd *m, + int i, int j, double coeff); + double MatrixXd_get_coeff (const struct C_MatrixXd *m, + int i, int j); + void MatrixXd_print (const struct C_MatrixXd *m); + void MatrixXd_add (const struct C_MatrixXd *m1, + const struct C_MatrixXd *m2, + struct C_MatrixXd *result); + void MatrixXd_multiply (const struct C_MatrixXd *m1, + const struct C_MatrixXd *m2, + struct C_MatrixXd *result); + + // the C_Map_MatrixXd class, wraps some of the functionality + // of Eigen::Map<MatrixXd> + struct C_Map_MatrixXd* Map_MatrixXd_new(double *array, int rows, int cols); + void Map_MatrixXd_delete (struct C_Map_MatrixXd *m); + void Map_MatrixXd_set_zero (struct C_Map_MatrixXd *m); + void Map_MatrixXd_copy (struct C_Map_MatrixXd *dst, + const struct C_Map_MatrixXd *src); + void Map_MatrixXd_copy_matrix(struct C_Map_MatrixXd *dst, + const struct C_MatrixXd *src); + void Map_MatrixXd_set_coeff (struct C_Map_MatrixXd *m, + int i, int j, double coeff); + double Map_MatrixXd_get_coeff (const struct C_Map_MatrixXd *m, + int i, int j); + void Map_MatrixXd_print (const struct C_Map_MatrixXd *m); + void Map_MatrixXd_add (const struct C_Map_MatrixXd *m1, + const struct C_Map_MatrixXd *m2, + struct C_Map_MatrixXd *result); + void Map_MatrixXd_multiply (const struct C_Map_MatrixXd *m1, + const struct C_Map_MatrixXd *m2, + struct C_Map_MatrixXd *result); + +#ifdef __cplusplus +} // end extern "C" +#endif \ No newline at end of file
diff --git a/demos/mix_eigen_and_c/example.c b/demos/mix_eigen_and_c/example.c new file mode 100644 index 0000000..4391a6d --- /dev/null +++ b/demos/mix_eigen_and_c/example.c
@@ -0,0 +1,80 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <http://www.gnu.org/licenses/>. + +#include "binary_library.h" +#include "stdio.h" + +void demo_MatrixXd() +{ + struct C_MatrixXd *matrix1, *matrix2, *result; + printf("*** demo_MatrixXd ***\n"); + + matrix1 = MatrixXd_new(3, 3); + MatrixXd_set_zero(matrix1); + MatrixXd_set_coeff(matrix1, 0, 1, 2.5); + MatrixXd_set_coeff(matrix1, 1, 0, 1.4); + printf("Here is matrix1:\n"); + MatrixXd_print(matrix1); + + matrix2 = MatrixXd_new(3, 3); + MatrixXd_multiply(matrix1, matrix1, matrix2); + printf("Here is matrix1*matrix1:\n"); + MatrixXd_print(matrix2); + + MatrixXd_delete(matrix1); + MatrixXd_delete(matrix2); +} + +// this helper function takes a plain C array and prints it in one line +void print_array(double *array, int n) +{ + struct C_Map_MatrixXd *m = Map_MatrixXd_new(array, 1, n); + Map_MatrixXd_print(m); + Map_MatrixXd_delete(m); +} + +void demo_Map_MatrixXd() +{ + struct C_Map_MatrixXd *map; + double array[5]; + int i; + printf("*** demo_Map_MatrixXd ***\n"); + + for(i = 0; i < 5; ++i) array[i] = i; + printf("Initially, the array is:\n"); + print_array(array, 5); + + map = Map_MatrixXd_new(array, 5, 1); + Map_MatrixXd_add(map, map, map); + Map_MatrixXd_delete(map); + + printf("Now the array is:\n"); + print_array(array, 5); +} + +int main() +{ + demo_MatrixXd(); + demo_Map_MatrixXd(); +}
diff --git a/demos/opengl/CMakeLists.txt b/demos/opengl/CMakeLists.txt index 968ed6c..b98a30c 100644 --- a/demos/opengl/CMakeLists.txt +++ b/demos/opengl/CMakeLists.txt
@@ -1,4 +1,3 @@ - find_package(Qt4 REQUIRED) find_package(OpenGL REQUIRED) @@ -14,6 +13,7 @@ qt4_automoc(${quaternion_demo_SRCS}) add_executable(quaternion_demo ${quaternion_demo_SRCS}) +add_dependencies(demos quaternion_demo) target_link_libraries(quaternion_demo ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY}
diff --git a/demos/opengl/camera.cpp b/demos/opengl/camera.cpp index a785caf..2659852 100644 --- a/demos/opengl/camera.cpp +++ b/demos/opengl/camera.cpp
@@ -260,7 +260,7 @@ Vector3f Camera::unProject(const Vector2f& uv, float depth) const { - Matrix4f inv = mViewMatrix.inverse(); + Matrix4f inv = mViewMatrix.inverse().matrix(); return unProject(uv, depth, inv); }
diff --git a/doc/C07_TutorialSparse.dox b/doc/C07_TutorialSparse.dox index 3a71828..ae96d77 100644 --- a/doc/C07_TutorialSparse.dox +++ b/doc/C07_TutorialSparse.dox
@@ -129,7 +129,7 @@ for (SparseVector<double>::InnerIterator it(vec); it; ++it) { it.value(); // == vec[ it.index() ] - it.index(); + it.index(); } \endcode </td></tr> @@ -138,58 +138,51 @@ \section TutorialSparseFilling Filling a sparse matrix -A DynamicSparseMatrix object can be set and updated just like any dense matrix using the coeffRef(row,col) method. If the coefficient is not stored yet, then it will be inserted in the matrix. Here is an example: +Owing to the special storage scheme of a SparseMatrix, it is obvious that for performance reasons a sparse matrix cannot be filled as easily as a dense matrix. For instance the cost of a purely random insertion into a SparseMatrix is in O(nnz) where nnz is the current number of non zeros. In order to cover all uses cases with best efficiency, Eigen provides various mechanisms, from the easiest but slowest, to the fastest but restrictive one. + +If you don't have any prior knowledge about the order your matrix will be filled, then the best choice is to use a DynamicSparseMatrix. With a DynamicSparseMatrix, you can add or modify any coefficients at any time using the coeffRef(row,col) method. Here is an example: \code DynamicSparseMatrix<float> aux(1000,1000); +aux.reserve(estimated_number_of_non_zero); // optional for (...) - for each i - for each j interacting with i - aux.coeffRef(i,j) += foo(o1,o2); -SparseMatrix<float> mat(aux); // convert the DynamicSparseMatrix to a SparseMatrix + for each j // the j can be random + for each i interacting with j // the i can be random + aux.coeffRef(i,j) += foo(i,j); +\endcode +Then the DynamicSparseMatrix object can be converted to a compact SparseMatrix to be used, e.g., by one of our supported solver: +\code +SparseMatrix<float> mat(aux); \endcode -Sometimes, however, we simply want to set all the coefficients of a matrix before using it through standard matrix operations (addition, product, etc.). In that case it faster to use the low-level startFill()/fill()/fillrand()/endFill() interface. Even though this interface is availabe for both sparse matrix types, their respective restrictions slightly differ from one representation to the other. In all case, a call to startFill() set the matrix to zero, and the fill*() functions will fail if the coefficient already exist. +In order to optimize this process, instead of the generic coeffRef(i,j) method one can also use: + - \code m.insert(i,j) = value; \endcode which assumes the coefficient of coordinate (row,col) does not already exist (otherwise this is a programming error and your program will stop). + - \code m.insertBack(i,j) = value; \endcode which, in addition to the requirements of insert(), also assumes that the coefficient of coordinate (row,col) will be inserted at the end of the target inner-vector. More precisely, if the matrix m is column major, then the row index of the last non zero coefficient of the j-th column must be smaller than i. -As a first difference, for SparseMatrix, the fill*() functions can only be called inside a startFill()/endFill() pair, and no other member functions are allowed during the filling process, i.e., until endFill() has been called. On the other hand, a DynamicSparseMatrix is always in a stable state, and the startFill()/endFill() functions are only for compatibility purpose. -Another difference is that the fill*() functions must be called with increasing outer indices for a SparseMatrix, while they can be random for a DynamicSparseMatrix. +Actually, the SparseMatrix class also supports random insertion via the insert() method. However, its uses should be reserved in cases where the inserted non zero is nearly the last one of the compact storage array. In practice, this means it should be used only to perform random (or sorted) insertion into the current inner-vector while filling the inner-vectors in an increasing order. Moreover, with a SparseMatrix an insertion session must be closed by a call to finalize() before any use of the matrix. Here is an example for a column major matrix: -Finally, the fill() function assumes the coefficient are inserted in a sorted order per inner vector, while the fillrand() variante allows random insertions (the outer indices must still be sorted for SparseMatrix). - -Some examples: - -1 - If you can set the coefficients in exactly the same order that the storage order, then the matrix can be filled directly and very efficiently. Here is an example initializing a random, row-major sparse matrix: \code -SparseMatrix<double,RowMajor> m(rows,cols); -m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional) -for (int i=0; i\<rows; ++i) - for (int j=0; j\<cols; ++j) - if (rand()\<percent_of_non_zero) - m.fill(i,j) = rand(); -m.endFill(); +SparseMatrix<float> mat(1000,1000); +mat.reserve(estimated_number_of_non_zero); // optional +for each j // should be in increasing order for performance reasons + for each i interacting with j // the i can be random + mat.insert(i,j) = foo(i,j); // optional for a DynamicSparseMatrix +mat.finalize(); \endcode -2 - If you can set each outer vector in a consistent order, but do not have sorted data for each inner vector, then you can use fillrand() instead of fill(): -\code -SparseMatrix<double,RowMajor> m(rows,cols); -m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional) -for (int i=0; i\<rows; ++i) - for (int k=0; k\<cols*percent_of_non_zero; ++k) - m.fillrand(i,rand(0,cols)) = rand(); -m.endFill(); -\endcode -The fillrand() function performs a sorted insertion into an array sequentially stored in memory and requires a copy of all coefficients stored after its target position. This method is therefore reserved for matrices having only a few elements per row/column (up to 50) and works better if the insertion are almost sorted. +Finally, the fastest way to fill a SparseMatrix object is to insert the elements in a purely coherence order (increasing inner index per increasing outer index). To this end, Eigen provides a very low but optimal API and illustrated below: -3 - Eventually, if none of the above solution is practicable for you, then you have to use a RandomSetter which temporarily wraps the matrix into a more flexible hash map allowing complete random accesses: \code -SparseMatrix<double,RowMajor> m(rows,cols); +SparseMatrix<float> mat(1000,1000); +mat.reserve(estimated_number_of_non_zero); // optional +for(int j=0; j<1000; ++j) { - RandomSetter<SparseMatrix<double,RowMajor> > setter(m); - for (int k=0; k\<cols*rows*percent_of_non_zero; ++k) - setter(rand(0,rows), rand(0,cols)) = rand(); + mat.startVec(j); // optional for a DynamicSparseMatrix + for each i interacting with j // with increasing i + mat.insertBack(i,j) = foo(i,j); } +mat.finalize(); // optional for a DynamicSparseMatrix \endcode -The matrix \c m is set at the destruction of the setter, hence the use of a nested block. This imposed syntax has the advantage to emphasize the critical section where m is not valid and cannot be used. \section TutorialSparseFeatureSet Supported operators and functions
diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 7738663..1d20a60 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt
@@ -1,3 +1,4 @@ +project(EigenDoc) set_directory_properties(PROPERTIES EXCLUDE_FROM_ALL TRUE)
diff --git a/doc/D11_UnalignedArrayAssert.dox b/doc/D11_UnalignedArrayAssert.dox index e95c767..226b204 100644 --- a/doc/D11_UnalignedArrayAssert.dox +++ b/doc/D11_UnalignedArrayAssert.dox
@@ -21,6 +21,7 @@ - \ref c3 - \ref c4 - \ref explanation + - \ref getrid \section where Where in my own code is the cause of the problem? @@ -101,6 +102,16 @@ However there are a few corner cases where these alignment settings get overridden: they are the possible causes for this assertion. +\section getrid I don't care about vectorization, how do I get rid of that stuff? + +Two possibilities: +<ul> + <li>Define EIGEN_DONT_ALIGN. That disables all 128-bit alignment code, and in particular everything vectorization-related. But do note that this in particular breaks ABI compatibility with vectorized code.</li> + <li>Or define both EIGEN_DONT_VECTORIZE and EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT. This keeps the 128-bit alignment code and thus preserves ABI compatibility.</li> +</ul> + +For more information, see <a href="http://eigen.tuxfamily.org/index.php?title=FAQ#I_disabled_vectorization.2C_but_I.27m_still_getting_annoyed_about_alignment_issues.21">this FAQ</a>. + */ }
diff --git a/doc/I00_CustomizingEigen.dox b/doc/I00_CustomizingEigen.dox index 8182a86..7c8ba9b 100644 --- a/doc/I00_CustomizingEigen.dox +++ b/doc/I00_CustomizingEigen.dox
@@ -125,7 +125,7 @@ 3 - define a couple of math functions for your type such as: ei_sqrt, ei_abs, etc... (see the file Eigen/src/Core/MathFunctions.h) -Here is a concrete example adding support for the Adolc's \c adouble type. <a href="http://www.math.tu-dresden.de/~adol-c/">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives. +Here is a concrete example adding support for the Adolc's \c adouble type. <a href="https://projects.coin-or.org/ADOL-C">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives. \code #ifndef ADLOCSUPPORT_H
diff --git a/scripts/eigen_gen_credits b/scripts/eigen_gen_credits new file mode 100755 index 0000000..c674167 --- /dev/null +++ b/scripts/eigen_gen_credits
@@ -0,0 +1,23 @@ +#!/bin/sh + +# this script must be run from the eigen2/ directory. +# when running hg churn from the scripts/ subdir, i hit a divide-by-zero error. +# +# like this: +# cd eigen2 +# USER=yourtuxfamilyuser scripts/eigen_gen_credits + +rm -f eigen_gen_credits.log + +hg pull >> eigen_gen_credits.log + +wget http://eigen.tuxfamily.org/index.php?title=ContributorsInfo -O online-info.out -a eigen_gen_credits.log +hg churn -r 37: --changesets -t {author} > churn-changesets.out +hg churn -r 37: -t {author} > churn-changedlines.out + +g++ scripts/eigen_gen_credits.cpp -o e + +./e > credits.out + +rsync credits.out $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/credits.out || (echo "upload failed"; exit 1) +ssh $USER@ssh.tuxfamily.org "cd eigen/eigen.tuxfamily.org-web/htdocs; chmod 664 credits.out; echo Main_Page | /usr/bin/php maintenance/purgeList.php"
diff --git a/scripts/eigen_gen_credits.cpp b/scripts/eigen_gen_credits.cpp new file mode 100644 index 0000000..086548e --- /dev/null +++ b/scripts/eigen_gen_credits.cpp
@@ -0,0 +1,213 @@ +#include <string> +#include <sstream> +#include <iostream> +#include <fstream> +#include <iomanip> +#include <map> +#include <list> + +using namespace std; + +// this function takes a line that may contain a name and/or email address, +// and returns just the name, while fixing the "bad cases". +std::string contributor_name(const std::string& line) +{ + string result; + size_t position_of_email_address = line.find_first_of('<'); + if(position_of_email_address != string::npos) + { + // there is an e-mail address. + + // Hauke once committed as "John Smith", fix that. + if(line.find("hauke.heibel") != string::npos) + result = "Hauke Heibel"; + else + { + // just remove the e-mail address + result = line.substr(0, position_of_email_address); + } + } + else + { + // there is no e-mail address. + + if(line.find("convert-repo") != string::npos) + result = ""; + else + result = line; + } + + // remove trailing spaces + size_t length = result.length(); + while(length >= 1 && result[length-1] == ' ') result.erase(--length); + + return result; +} + +// parses hg churn output to generate a contributors map. +map<string,int> contributors_map_from_churn_output(const char *filename) +{ + map<string,int> contributors_map; + + string line; + ifstream churn_out; + churn_out.open(filename, ios::in); + while(!getline(churn_out,line).eof()) + { + // remove the histograms "******" that hg churn may draw at the end of some lines + size_t first_star = line.find_first_of('*'); + if(first_star != string::npos) line.erase(first_star); + + // remove trailing spaces + size_t length = line.length(); + while(length >= 1 && line[length-1] == ' ') line.erase(--length); + + // now the last space indicates where the number starts + size_t last_space = line.find_last_of(' '); + + // get the number (of changesets or of modified lines for each contributor) + int number; + istringstream(line.substr(last_space+1)) >> number; + + // get the name of the contributor + line.erase(last_space); + string name = contributor_name(line); + + map<string,int>::iterator it = contributors_map.find(name); + // if new contributor, insert + if(it == contributors_map.end()) + contributors_map.insert(pair<string,int>(name, number)); + // if duplicate, just add the number + else + it->second += number; + } + churn_out.close(); + + return contributors_map; +} + +// find the last name, i.e. the last word. +// for "van den Schbling" types of last names, that's not a problem, that's actually what we want. +string lastname(const string& name) +{ + size_t last_space = name.find_last_of(' '); + if(last_space >= name.length()-1) return name; + else return name.substr(last_space+1); +} + +struct contributor +{ + string name; + int changedlines; + int changesets; + string url; + string misc; + + contributor() : changedlines(0), changesets(0) {} + + bool operator < (const contributor& other) + { + return lastname(name).compare(lastname(other.name)) < 0; + } +}; + +void add_online_info_into_contributors_list(list<contributor>& contributors_list, const char *filename) +{ + string line; + ifstream online_info; + online_info.open(filename, ios::in); + while(!getline(online_info,line).eof()) + { + string hgname, realname, url, misc; + + size_t last_bar = line.find_last_of('|'); + if(last_bar == string::npos) continue; + if(last_bar < line.length()) + misc = line.substr(last_bar+1); + line.erase(last_bar); + + last_bar = line.find_last_of('|'); + if(last_bar == string::npos) continue; + if(last_bar < line.length()) + url = line.substr(last_bar+1); + line.erase(last_bar); + + last_bar = line.find_last_of('|'); + if(last_bar == string::npos) continue; + if(last_bar < line.length()) + realname = line.substr(last_bar+1); + line.erase(last_bar); + + hgname = line; + + // remove the example line + if(hgname.find("MercurialName") != string::npos) continue; + + list<contributor>::iterator it; + for(it=contributors_list.begin(); it != contributors_list.end() && it->name != hgname; ++it) + {} + + if(it == contributors_list.end()) + { + contributor c; + c.name = realname; + c.url = url; + c.misc = misc; + contributors_list.push_back(c); + } + else + { + it->name = realname; + it->url = url; + it->misc = misc; + } + } +} + +int main() +{ + // parse the hg churn output files + map<string,int> contributors_map_for_changedlines = contributors_map_from_churn_output("churn-changedlines.out"); + //map<string,int> contributors_map_for_changesets = contributors_map_from_churn_output("churn-changesets.out"); + + // merge into the contributors list + list<contributor> contributors_list; + map<string,int>::iterator it; + for(it=contributors_map_for_changedlines.begin(); it != contributors_map_for_changedlines.end(); ++it) + { + contributor c; + c.name = it->first; + c.changedlines = it->second; + c.changesets = 0; //contributors_map_for_changesets.find(it->first)->second; + contributors_list.push_back(c); + } + + add_online_info_into_contributors_list(contributors_list, "online-info.out"); + + contributors_list.sort(); + + cout << "{| cellpadding=\"5\"\n"; + cout << "!\n"; + cout << "! Lines changed\n"; + cout << "!\n"; + + list<contributor>::iterator itc; + int i = 0; + for(itc=contributors_list.begin(); itc != contributors_list.end(); ++itc) + { + if(itc->name.length() == 0) continue; + if(i%2) cout << "|-\n"; + else cout << "|- style=\"background:#FFFFD0\"\n"; + if(itc->url.length()) + cout << "| [" << itc->url << " " << itc->name << "]\n"; + else + cout << "| " << itc->name << "\n"; + if(itc->changedlines) + cout << "| " << itc->changedlines << "\n"; + else + cout << "| (no information)\n"; + cout << "| " << itc->misc << "\n"; + i++; + } + cout << "|}" << endl; +}
diff --git a/scripts/eigen_gen_docs b/scripts/eigen_gen_docs index 3cdacc1..17f5736 100644 --- a/scripts/eigen_gen_docs +++ b/scripts/eigen_gen_docs
@@ -1,18 +1,17 @@ #!/bin/sh # configuration -USER='orzel' +# You should call this script with USER set as you want, else some default +# will be used +USER=${USER:-'orzel'} -# step 1 : update -hg pull -u || (echo "update failed"; exit 1) - -# step 2 : build +# step 1 : build # todo if 'build is not there, create one: #mkdir build (cd build && cmake .. && make -j3 doc) || (echo "make failed"; exit 1) #todo: n+1 where n = number of cpus -#step 3 : upload +#step 2 : upload BRANCH=`hg branch` if [ $BRANCH == "default" ] then
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4e279ea..b469714 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt
@@ -1,6 +1,6 @@ - +project(EigenTesting) +add_custom_target(btest) include(EigenTesting) -enable_testing() ei_init_testing() find_package(GSL) @@ -94,6 +94,7 @@ ei_add_test(linearstructure) ei_add_test(cwiseop) ei_add_test(redux) +ei_add_test(visitor) ei_add_test(product_small) ei_add_test(product_large ${EI_OFLAG}) ei_add_test(product_extra ${EI_OFLAG})
diff --git a/test/geo_homogeneous.cpp b/test/geo_homogeneous.cpp index c4db32f..4b66e48 100644 --- a/test/geo_homogeneous.cpp +++ b/test/geo_homogeneous.cpp
@@ -97,7 +97,7 @@ pts.setRandom(Size,5); Rt_pts1 = Rt * pts.colwise().homogeneous(); - std::cerr << (Rt_pts1 - pts).sum() << "\n"; + // std::cerr << (Rt_pts1 - pts).sum() << "\n"; VERIFY_IS_MUCH_SMALLER_THAN( (Rt_pts1 - pts).sum(), Scalar(1)); }
diff --git a/test/geo_transformations.cpp b/test/geo_transformations.cpp index 914dbaf..23b2973 100644 --- a/test/geo_transformations.cpp +++ b/test/geo_transformations.cpp
@@ -102,7 +102,14 @@ a = ei_random<Scalar>(-Scalar(0.4)*Scalar(M_PI), Scalar(0.4)*Scalar(M_PI)); q1 = AngleAxisx(a, v0.normalized()); Transform3 t0, t1, t2; + + // first test setIdentity() and Identity() t0.setIdentity(); + VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity()); + t0.matrix().setZero(); + t0 = Transform3::Identity(); + VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity()); + t0.linear() = q1.toRotationMatrix(); t1.setIdentity(); t1.linear() = q1.toRotationMatrix(); @@ -296,10 +303,10 @@ t0.setIdentity(); t0.translate(v0); t0.linear().setRandom(); - VERIFY_IS_APPROX(t0.inverse(Affine), t0.matrix().inverse()); + VERIFY_IS_APPROX(t0.inverse(Affine).matrix(), t0.matrix().inverse()); t0.setIdentity(); t0.translate(v0).rotate(q1); - VERIFY_IS_APPROX(t0.inverse(Isometry), t0.matrix().inverse()); + VERIFY_IS_APPROX(t0.inverse(Isometry).matrix(), t0.matrix().inverse()); } // test extract rotation and aligned scaling
diff --git a/test/householder.cpp b/test/householder.cpp index 7d30089..b272794 100644 --- a/test/householder.cpp +++ b/test/householder.cpp
@@ -43,7 +43,7 @@ Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols)); Scalar* tmp = &_tmp.coeffRef(0,0); - + Scalar beta; RealScalar alpha; EssentialVectorType essential; @@ -58,7 +58,7 @@ v2 = v1; v1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(v1.norm(), v2.norm()); - + MatrixType m1(rows, cols), m2(rows, cols); @@ -72,7 +72,7 @@ VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm()); VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m1(0,0)), ei_real(m1(0,0))); VERIFY_IS_APPROX(ei_real(m1(0,0)), alpha); - + v1 = VectorType::Random(rows); if(even) v1.end(rows-1).setZero(); SquareMatrixType m3(rows,rows), m4(rows,rows); @@ -84,6 +84,9 @@ VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm()); VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0))); VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha); + + // test householder sequence + // TODO test HouseholderSequence } void test_householder()
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 5940b84..2e3f089 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp
@@ -36,14 +36,14 @@ RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime }; - + typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType; - + MatrixType a; if(pickrandom) a = MatrixType::Random(rows,cols); else a = m; @@ -53,7 +53,7 @@ sigma.diagonal() = svd.singularValues().template cast<Scalar>(); MatrixUType u = svd.matrixU(); MatrixVType v = svd.matrixV(); - + VERIFY_IS_APPROX(a, u * sigma * v.adjoint()); VERIFY_IS_UNITARY(u); VERIFY_IS_UNITARY(v); @@ -98,7 +98,7 @@ } CALL_SUBTEST(( svd<MatrixXf,0>(MatrixXf(300,200)) )); CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyColsAsRows>(MatrixXcd(100,150)) )); - + CALL_SUBTEST(( svd_verify_assert<Matrix3f>() )); CALL_SUBTEST(( svd_verify_assert<Matrix3d>() )); CALL_SUBTEST(( svd_verify_assert<MatrixXf>() ));
diff --git a/test/main.h b/test/main.h index 619fc9e..51b7198 100644 --- a/test/main.h +++ b/test/main.h
@@ -40,6 +40,11 @@ #define DEFAULT_REPEAT 10 +#ifdef __ICC +// disable warning #279: controlling expression is constant +#pragma warning disable 279 +#endif + namespace Eigen { static std::vector<std::string> g_test_stack; @@ -243,18 +248,22 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m) { typedef typename ei_traits<MatrixType>::Scalar Scalar; - typedef Matrix<Scalar, Dynamic, 1> VectorType; + enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; - MatrixType a = MatrixType::Random(rows,rows); + typedef Matrix<Scalar, Dynamic, 1> VectorType; + typedef Matrix<Scalar, Rows, Rows> MatrixAType; + typedef Matrix<Scalar, Cols, Cols> MatrixBType; + + MatrixAType a = MatrixAType::Random(rows,rows); MatrixType d = MatrixType::Identity(rows,cols); - MatrixType b = MatrixType::Random(cols,cols); + MatrixBType b = MatrixBType::Random(cols,cols); // set the diagonal such that only desired_rank non-zero entries reamain const int diag_size = std::min(d.rows(),d.cols()); d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank); - HouseholderQR<MatrixType> qra(a); - HouseholderQR<MatrixType> qrb(b); + HouseholderQR<MatrixAType> qra(a); + HouseholderQR<MatrixBType> qrb(b); m = qra.matrixQ() * d * qrb.matrixQ(); }
diff --git a/test/packetmath.cpp b/test/packetmath.cpp index d86d40d..1745ae5 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp
@@ -99,10 +99,10 @@ const int PacketSize = ei_packet_traits<Scalar>::size; const int size = PacketSize*4; - EIGEN_ALIGN_128 Scalar data1[ei_packet_traits<Scalar>::size*4]; - EIGEN_ALIGN_128 Scalar data2[ei_packet_traits<Scalar>::size*4]; - EIGEN_ALIGN_128 Packet packets[PacketSize*2]; - EIGEN_ALIGN_128 Scalar ref[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Scalar data1[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Scalar data2[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Packet packets[PacketSize*2]; + EIGEN_ALIGN16 Scalar ref[ei_packet_traits<Scalar>::size*4]; for (int i=0; i<size; ++i) { data1[i] = ei_random<Scalar>(); @@ -202,9 +202,9 @@ const int PacketSize = ei_packet_traits<Scalar>::size; const int size = PacketSize*4; - EIGEN_ALIGN_128 Scalar data1[ei_packet_traits<Scalar>::size*4]; - EIGEN_ALIGN_128 Scalar data2[ei_packet_traits<Scalar>::size*4]; - EIGEN_ALIGN_128 Scalar ref[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Scalar data1[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Scalar data2[ei_packet_traits<Scalar>::size*4]; + EIGEN_ALIGN16 Scalar ref[ei_packet_traits<Scalar>::size*4]; for (int i=0; i<size; ++i) {
diff --git a/test/product_extra.cpp b/test/product_extra.cpp index 3ad99fc..8e55c60 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp
@@ -114,7 +114,6 @@ VERIFY_IS_APPROX(m1.col(j2).adjoint() * m1.block(0,j,m1.rows(),c), m1.col(j2).adjoint().eval() * m1.block(0,j,m1.rows(),c).eval()); VERIFY_IS_APPROX(m1.block(i,0,r,m1.cols()) * m1.row(i2).adjoint(), m1.block(i,0,r,m1.cols()).eval() * m1.row(i2).adjoint().eval()); - } void test_product_extra()
diff --git a/test/product_small.cpp b/test/product_small.cpp index 3aed5cf..182af71 100644 --- a/test/product_small.cpp +++ b/test/product_small.cpp
@@ -34,4 +34,10 @@ CALL_SUBTEST( product(Matrix4d()) ); CALL_SUBTEST( product(Matrix4f()) ); } + + { + // test compilation of (outer_product) * vector + Vector3f v = Vector3f::Random(); + VERIFY_IS_APPROX( (v * v.transpose()) * v, (v * v.transpose()).eval() * v); + } }
diff --git a/test/qr.cpp b/test/qr.cpp index f2e2eda..036a3c9 100644 --- a/test/qr.cpp +++ b/test/qr.cpp
@@ -41,20 +41,26 @@ for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); VERIFY_IS_APPROX(a, qrOfA.matrixQ() * r); +} - SquareMatrixType b = a.adjoint() * a; +template<typename MatrixType, int Cols2> void qr_fixedsize() +{ + enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; + typedef typename MatrixType::Scalar Scalar; + Matrix<Scalar,Rows,Cols> m1 = Matrix<Scalar,Rows,Cols>::Random(); + HouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1); - // check tridiagonalization - Tridiagonalization<SquareMatrixType> tridiag(b); - VERIFY_IS_APPROX(b, tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint()); + Matrix<Scalar,Rows,Cols> r = qr.matrixQR(); + // FIXME need better way to construct trapezoid + for(int i = 0; i < Rows; i++) for(int j = 0; j < Cols; j++) if(i>j) r(i,j) = Scalar(0); - // check hessenberg decomposition - HessenbergDecomposition<SquareMatrixType> hess(b); - VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); - VERIFY_IS_APPROX(tridiag.matrixT(), hess.matrixH()); - b = SquareMatrixType::Random(cols,cols); - hess.compute(b); - VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); + VERIFY_IS_APPROX(m1, qr.matrixQ() * r); + + Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2); + Matrix<Scalar,Rows,Cols2> m3 = m1*m2; + m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2); + qr.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); } template<typename MatrixType> void qr_invertible() @@ -78,7 +84,7 @@ m3 = MatrixType::Random(size,size); qr.solve(m3, &m2); VERIFY_IS_APPROX(m3, m1*m2); - + // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>(); @@ -105,11 +111,11 @@ void test_qr() { for(int i = 0; i < 1; i++) { - // FIXME : very weird bug here -// CALL_SUBTEST( qr(Matrix2f()) ); - CALL_SUBTEST( qr(Matrix4d()) ); - CALL_SUBTEST( qr(MatrixXf(47,40)) ); - CALL_SUBTEST( qr(MatrixXcd(17,7)) ); + CALL_SUBTEST( qr(MatrixXf(47,40)) ); + CALL_SUBTEST( qr(MatrixXcd(17,7)) ); + CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,4>, 2 >() )); + CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 4 >() )); + CALL_SUBTEST(( qr_fixedsize<Matrix<double,2,5>, 7 >() )); } for(int i = 0; i < g_repeat; i++) {
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 2838554..4b6f7dd 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp
@@ -42,18 +42,18 @@ VERIFY(!qr.isInjective()); VERIFY(!qr.isInvertible()); VERIFY(!qr.isSurjective()); - + MatrixType r = qr.matrixQR(); // FIXME need better way to construct trapezoid for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); - + MatrixType b = qr.matrixQ() * r; MatrixType c = MatrixType::Zero(rows,cols); - + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); VERIFY_IS_APPROX(m1, c); - + MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); @@ -63,6 +63,40 @@ VERIFY(!qr.solve(m3, &m2)); } +template<typename MatrixType, int Cols2> void qr_fixedsize() +{ + enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; + typedef typename MatrixType::Scalar Scalar; + int rank = ei_random<int>(1, std::min(int(Rows), int(Cols))-1); + Matrix<Scalar,Rows,Cols> m1; + createRandomMatrixOfRank(rank,Rows,Cols,m1); + ColPivotingHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1); + VERIFY_IS_APPROX(rank, qr.rank()); + VERIFY(Cols - qr.rank() == qr.dimensionOfKernel()); + VERIFY(!qr.isInjective()); + VERIFY(!qr.isInvertible()); + VERIFY(!qr.isSurjective()); + + Matrix<Scalar,Rows,Cols> r = qr.matrixQR(); + // FIXME need better way to construct trapezoid + for(int i = 0; i < Rows; i++) for(int j = 0; j < Cols; j++) if(i>j) r(i,j) = Scalar(0); + + Matrix<Scalar,Rows,Cols> b = qr.matrixQ() * r; + + Matrix<Scalar,Rows,Cols> c = MatrixType::Zero(Rows,Cols); + + for(int i = 0; i < Cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); + VERIFY_IS_APPROX(m1, c); + + Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2); + Matrix<Scalar,Rows,Cols2> m3 = m1*m2; + m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2); + VERIFY(qr.solve(m3, &m2)); + VERIFY_IS_APPROX(m3, m1*m2); + m3 = Matrix<Scalar,Rows,Cols2>::Random(Rows,Cols2); + VERIFY(!qr.solve(m3, &m2)); +} + template<typename MatrixType> void qr_invertible() { typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; @@ -116,12 +150,12 @@ void test_qr_colpivoting() { - for(int i = 0; i < 1; i++) { - // FIXME : very weird bug here -// CALL_SUBTEST( qr(Matrix2f()) ); + for(int i = 0; i < 1; i++) { CALL_SUBTEST( qr<MatrixXf>() ); CALL_SUBTEST( qr<MatrixXd>() ); CALL_SUBTEST( qr<MatrixXcd>() ); + CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,5>, 4 >() )); + CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 3 >() )); } for(int i = 0; i < g_repeat; i++) {
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index 525c669..3a37bcb 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp
@@ -46,14 +46,14 @@ MatrixType r = qr.matrixQR(); // FIXME need better way to construct trapezoid for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); - + MatrixType b = qr.matrixQ() * r; MatrixType c = MatrixType::Zero(rows,cols); - + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); VERIFY_IS_APPROX(m1, c); - + MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); @@ -88,7 +88,7 @@ m3 = MatrixType::Random(size,size); VERIFY(qr.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); - + // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
diff --git a/test/redux.cpp b/test/redux.cpp index 2a0dc97..951b34b 100644 --- a/test/redux.cpp +++ b/test/redux.cpp
@@ -120,7 +120,7 @@ CALL_SUBTEST( matrixRedux(MatrixXi(8, 12)) ); } for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( vectorRedux(VectorXf(5)) ); + CALL_SUBTEST( vectorRedux(Vector4f()) ); CALL_SUBTEST( vectorRedux(VectorXd(10)) ); CALL_SUBTEST( vectorRedux(VectorXf(33)) ); }
diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index b8fbf52..ed72bb7 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp
@@ -33,10 +33,24 @@ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; + // Check the basic machine-dependent constants. + { + int ibeta, it, iemin, iemax; + + ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers + it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa + iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent + iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent + + VERIFY( (!(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) || (it<=4 && ibeta <= 3 ) || it<2)) + && "the stable norm algorithm cannot be guaranteed on this computer"); + } + + int rows = m.rows(); int cols = m.cols(); - Scalar big = ei_random<Scalar>() * std::numeric_limits<RealScalar>::max() * RealScalar(1e-4); + Scalar big = ei_abs(ei_random<Scalar>()) * (std::numeric_limits<RealScalar>::max() * RealScalar(1e-4)); Scalar small = static_cast<RealScalar>(1)/big; MatrixType vzero = MatrixType::Zero(rows, cols),
diff --git a/test/testsuite.cmake b/test/testsuite.cmake index 37ee875..5d0cb65 100644 --- a/test/testsuite.cmake +++ b/test/testsuite.cmake
@@ -148,7 +148,7 @@ # which ctest command to use for running the dashboard SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}") # what cmake command to use for configuring this dashboard -SET (CTEST_CMAKE_COMMAND "${EIGEN_CMAKE_DIR}cmake -DEIGEN_BUILD_TESTS=on ") +SET (CTEST_CMAKE_COMMAND "${EIGEN_CMAKE_DIR}cmake -DEIGEN_CMAKE_RUN_FROM_CTEST=ON") #################################################################### # The values in this section are optional you can either @@ -175,9 +175,9 @@ else(EIGEN_GENERATOR_TYPE) set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake") SET (CTEST_INITIAL_CACHE " - MAKECOMMAND:STRING=nmake -i + MAKECOMMAND:STRING=nmake /i CMAKE_MAKE_PROGRAM:FILEPATH=nmake - CMAKE_GENERATOR:INTERNAL=NMake Makefiles + CMAKE_GENERATOR:INTERNAL=NMake Makefiles CMAKE_BUILD_TYPE:STRING=Release BUILDNAME:STRING=${EIGEN_BUILD_STRING} SITE:STRING=${EIGEN_SITE}
diff --git a/test/unalignedassert.cpp b/test/unalignedassert.cpp index ade1ab2..233268d 100644 --- a/test/unalignedassert.cpp +++ b/test/unalignedassert.cpp
@@ -24,52 +24,38 @@ #include "main.h" -struct Good1 +struct TestNew1 { MatrixXd m; // good: m will allocate its own array, taking care of alignment. - Good1() : m(20,20) {} + TestNew1() : m(20,20) {} }; -struct Good2 +struct TestNew2 { - Matrix3d m; // good: m's size isn't a multiple of 16 bytes, so m doesn't have to be aligned + Matrix3d m; // good: m's size isn't a multiple of 16 bytes, so m doesn't have to be 16-byte aligned, + // 8-byte alignment is good enough here, which we'll get automatically }; -struct Good3 +struct TestNew3 { - Vector2f m; // good: same reason + Vector2f m; // good: m's size isn't a multiple of 16 bytes, so m doesn't have to be 16-byte aligned }; -struct Bad4 -{ - Vector2d m; // bad: sizeof(m)%16==0 so alignment is required -}; - -struct Bad5 -{ - Matrix<float, 2, 6> m; // bad: same reason -}; - -struct Bad6 -{ - Matrix<double, 3, 4> m; // bad: same reason -}; - -struct Good7 +struct TestNew4 { EIGEN_MAKE_ALIGNED_OPERATOR_NEW Vector2d m; float f; // make the struct have sizeof%16!=0 to make it a little more tricky when we allow an array of 2 such objects }; -struct Good8 +struct TestNew5 { EIGEN_MAKE_ALIGNED_OPERATOR_NEW - float f; // try the f at first -- the EIGEN_ALIGN_128 attribute of m should make that still work + float f; // try the f at first -- the EIGEN_ALIGN16 attribute of m should make that still work Matrix4f m; }; -struct Good9 +struct TestNew6 { Matrix<float,2,2,DontAlign> m; // good: no alignment requested float f; @@ -94,34 +80,56 @@ #if EIGEN_ALIGN template<typename T> -void check_unalignedassert_bad() +void construct_at_boundary(int boundary) { - float buf[sizeof(T)+16]; - float *unaligned = buf; - while((reinterpret_cast<size_t>(unaligned)&0xf)==0) ++unaligned; // make sure unaligned is really unaligned - T *x = ::new(static_cast<void*>(unaligned)) T; + char buf[sizeof(T)+256]; + size_t _buf = reinterpret_cast<size_t>(buf); + _buf += (16 - (_buf % 16)); // make 16-byte aligned + _buf += boundary; // make exact boundary-aligned + T *x = ::new(reinterpret_cast<void*>(_buf)) T; x->~T(); } #endif void unalignedassert() { - check_unalignedassert_good<Good1>(); - check_unalignedassert_good<Good2>(); - check_unalignedassert_good<Good3>(); -#if EIGEN_ALIGN - VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad4>()); - VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad5>()); - VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Bad6>()); -#endif + construct_at_boundary<Vector2f>(4); + construct_at_boundary<Vector3f>(4); + construct_at_boundary<Vector4f>(16); + construct_at_boundary<Matrix2f>(16); + construct_at_boundary<Matrix3f>(4); + construct_at_boundary<Matrix4f>(16); + + construct_at_boundary<Vector2d>(16); + construct_at_boundary<Vector3d>(4); + construct_at_boundary<Vector4d>(16); + construct_at_boundary<Matrix2d>(16); + construct_at_boundary<Matrix3d>(4); + construct_at_boundary<Matrix4d>(16); + + construct_at_boundary<Vector2cf>(16); + construct_at_boundary<Vector3cf>(4); + construct_at_boundary<Vector2cd>(16); + construct_at_boundary<Vector3cd>(16); + + check_unalignedassert_good<TestNew1>(); + check_unalignedassert_good<TestNew2>(); + check_unalignedassert_good<TestNew3>(); - check_unalignedassert_good<Good7>(); - check_unalignedassert_good<Good8>(); - check_unalignedassert_good<Good9>(); + check_unalignedassert_good<TestNew4>(); + check_unalignedassert_good<TestNew5>(); + check_unalignedassert_good<TestNew6>(); check_unalignedassert_good<Depends<true> >(); - + #if EIGEN_ALIGN - VERIFY_RAISES_ASSERT(check_unalignedassert_bad<Depends<false> >()); + VERIFY_RAISES_ASSERT(construct_at_boundary<Vector4f>(8)); + VERIFY_RAISES_ASSERT(construct_at_boundary<Matrix4f>(8)); + VERIFY_RAISES_ASSERT(construct_at_boundary<Vector2d>(8)); + VERIFY_RAISES_ASSERT(construct_at_boundary<Vector4d>(8)); + VERIFY_RAISES_ASSERT(construct_at_boundary<Matrix2d>(8)); + VERIFY_RAISES_ASSERT(construct_at_boundary<Matrix4d>(8)); + VERIFY_RAISES_ASSERT(construct_at_boundary<Vector2cf>(8)); + VERIFY_RAISES_ASSERT(construct_at_boundary<Vector2cd>(8)); #endif }
diff --git a/test/visitor.cpp b/test/visitor.cpp new file mode 100644 index 0000000..b78782b --- /dev/null +++ b/test/visitor.cpp
@@ -0,0 +1,131 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// 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 <http://www.gnu.org/licenses/>. + +#include "main.h" + +template<typename MatrixType> void matrixVisitor(const MatrixType& p) +{ + typedef typename MatrixType::Scalar Scalar; + + int rows = p.rows(); + int cols = p.cols(); + + // construct a random matrix where all coefficients are different + MatrixType m; + m = MatrixType::Random(rows, cols); + for(int i = 0; i < m.size(); i++) + for(int i2 = 0; i2 < i; i2++) + while(m(i) == m(i2)) // yes, == + m(i) = ei_random<Scalar>(); + + Scalar minc = Scalar(1000), maxc = Scalar(-1000); + int minrow,mincol,maxrow,maxcol; + for(int j = 0; j < cols; j++) + for(int i = 0; i < rows; i++) + { + if(m(i,j) < minc) + { + minc = m(i,j); + minrow = i; + mincol = j; + } + if(m(i,j) > maxc) + { + maxc = m(i,j); + maxrow = i; + maxcol = j; + } + } + int eigen_minrow, eigen_mincol, eigen_maxrow, eigen_maxcol; + Scalar eigen_minc, eigen_maxc; + eigen_minc = m.minCoeff(&eigen_minrow,&eigen_mincol); + eigen_maxc = m.maxCoeff(&eigen_maxrow,&eigen_maxcol); + VERIFY(minrow == eigen_minrow); + VERIFY(maxrow == eigen_maxrow); + VERIFY(mincol == eigen_mincol); + VERIFY(maxcol == eigen_maxcol); + VERIFY_IS_APPROX(minc, eigen_minc); + VERIFY_IS_APPROX(maxc, eigen_maxc); + VERIFY_IS_APPROX(minc, m.minCoeff()); + VERIFY_IS_APPROX(maxc, m.maxCoeff()); +} + +template<typename VectorType> void vectorVisitor(const VectorType& w) +{ + typedef typename VectorType::Scalar Scalar; + + int size = w.size(); + + // construct a random vector where all coefficients are different + VectorType v; + v = VectorType::Random(size); + for(int i = 0; i < size; i++) + for(int i2 = 0; i2 < i; i2++) + while(v(i) == v(i2)) // yes, == + v(i) = ei_random<Scalar>(); + + Scalar minc = Scalar(1000), maxc = Scalar(-1000); + int minidx,maxidx; + for(int i = 0; i < size; i++) + { + if(v(i) < minc) + { + minc = v(i); + minidx = i; + } + if(v(i) > maxc) + { + maxc = v(i); + maxidx = i; + } + } + int eigen_minidx, eigen_maxidx; + Scalar eigen_minc, eigen_maxc; + eigen_minc = v.minCoeff(&eigen_minidx); + eigen_maxc = v.maxCoeff(&eigen_maxidx); + VERIFY(minidx == eigen_minidx); + VERIFY(maxidx == eigen_maxidx); + VERIFY_IS_APPROX(minc, eigen_minc); + VERIFY_IS_APPROX(maxc, eigen_maxc); + VERIFY_IS_APPROX(minc, v.minCoeff()); + VERIFY_IS_APPROX(maxc, v.maxCoeff()); +} + +void test_visitor() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( matrixVisitor(Matrix<float, 1, 1>()) ); + CALL_SUBTEST( matrixVisitor(Matrix2f()) ); + CALL_SUBTEST( matrixVisitor(Matrix4d()) ); + CALL_SUBTEST( matrixVisitor(MatrixXd(8, 12)) ); + CALL_SUBTEST( matrixVisitor(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) ); + CALL_SUBTEST( matrixVisitor(MatrixXi(8, 12)) ); + } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( vectorVisitor(Vector4f()) ); + CALL_SUBTEST( vectorVisitor(VectorXd(10)) ); + CALL_SUBTEST( vectorVisitor(RowVectorXd(10)) ); + CALL_SUBTEST( vectorVisitor(VectorXf(33)) ); + } +}
diff --git a/unsupported/CMakeLists.txt b/unsupported/CMakeLists.txt index 895fcdb..7104085 100644 --- a/unsupported/CMakeLists.txt +++ b/unsupported/CMakeLists.txt
@@ -1,9 +1,3 @@ - add_subdirectory(Eigen) - -add_subdirectory(doc) - -if(EIGEN_BUILD_TESTS) - add_subdirectory(test) -endif(EIGEN_BUILD_TESTS) - +add_subdirectory(doc EXCLUDE_FROM_ALL) +add_subdirectory(test EXCLUDE_FROM_ALL)
diff --git a/unsupported/Eigen/AdolcForward b/unsupported/Eigen/AdolcForward index 9a54a3a..70aa478 100644 --- a/unsupported/Eigen/AdolcForward +++ b/unsupported/Eigen/AdolcForward
@@ -29,7 +29,7 @@ // // This file provides support for adolc's adouble type in forward mode. // ADOL-C is a C++ automatic differentiation library, -// see http://www.math.tu-dresden.de/~adol-c/ for more information. +// see https://projects.coin-or.org/ADOL-C for more information. // // Note that the maximal number of directions is controlled by // the preprocessor token NUMBER_DIRECTIONS. The default is 2. @@ -63,7 +63,7 @@ * \defgroup AdolcForward_Module Adolc forward module * This module provides support for adolc's adouble type in forward mode. * ADOL-C is a C++ automatic differentiation library, - * see http://www.math.tu-dresden.de/~adol-c/ for more information. + * see https://projects.coin-or.org/ADOL-C for more information. * It mainly consists in: * - a struct Eigen::NumTraits<adtl::adouble> specialization * - overloads of ei_* math function for adtl::adouble type.
diff --git a/unsupported/Eigen/AlignedVector3 b/unsupported/Eigen/AlignedVector3 index aaec5f9..f20fad6 100644 --- a/unsupported/Eigen/AlignedVector3 +++ b/unsupported/Eigen/AlignedVector3
@@ -63,37 +63,37 @@ typedef Matrix<_Scalar,4,1> CoeffType; CoeffType m_coeffs; public: - + EIGEN_GENERIC_PUBLIC_INTERFACE(AlignedVector3) using Base::operator*; - + inline int rows() const { return 3; } inline int cols() const { return 1; } - + inline const Scalar& coeff(int row, int col) const { return m_coeffs.coeff(row, col); } - + inline Scalar& coeffRef(int row, int col) { return m_coeffs.coeffRef(row, col); } - + inline const Scalar& coeff(int index) const { return m_coeffs.coeff(index); } inline Scalar& coeffRef(int index) { return m_coeffs.coeffRef(index);} - - + + inline AlignedVector3(const Scalar& x, const Scalar& y, const Scalar& z) : m_coeffs(x, y, z, Scalar(0)) {} - + inline AlignedVector3(const AlignedVector3& other) - : m_coeffs(other.m_coeffs) + : Base(), m_coeffs(other.m_coeffs) {} - + template<typename XprType, int Size=XprType::SizeAtCompileTime> struct generic_assign_selector {}; - + template<typename XprType> struct generic_assign_selector<XprType,4> { inline static void run(AlignedVector3& dest, const XprType& src) @@ -101,7 +101,7 @@ dest.m_coeffs = src; } }; - + template<typename XprType> struct generic_assign_selector<XprType,3> { inline static void run(AlignedVector3& dest, const XprType& src) @@ -110,44 +110,44 @@ dest.m_coeffs.w() = Scalar(0); } }; - + template<typename Derived> inline explicit AlignedVector3(const MatrixBase<Derived>& other) { generic_assign_selector<Derived>::run(*this,other.derived()); } - + inline AlignedVector3& operator=(const AlignedVector3& other) { m_coeffs = other.m_coeffs; return *this; } - - + + inline AlignedVector3 operator+(const AlignedVector3& other) const { return AlignedVector3(m_coeffs + other.m_coeffs); } - + inline AlignedVector3& operator+=(const AlignedVector3& other) { m_coeffs += other.m_coeffs; return *this; } - + inline AlignedVector3 operator-(const AlignedVector3& other) const { return AlignedVector3(m_coeffs - other.m_coeffs); } - + inline AlignedVector3 operator-=(const AlignedVector3& other) { m_coeffs -= other.m_coeffs; return *this; } - + inline AlignedVector3 operator*(const Scalar& s) const { return AlignedVector3(m_coeffs * s); } - + inline friend AlignedVector3 operator*(const Scalar& s,const AlignedVector3& vec) { return AlignedVector3(s * vec.m_coeffs); } - + inline AlignedVector3& operator*=(const Scalar& s) { m_coeffs *= s; return *this; } - + inline AlignedVector3 operator/(const Scalar& s) const { return AlignedVector3(m_coeffs / s); } - + inline AlignedVector3& operator/=(const Scalar& s) { m_coeffs /= s; return *this; } - + inline Scalar dot(const AlignedVector3& other) const { ei_assert(m_coeffs.w()==Scalar(0)); @@ -164,29 +164,29 @@ { return AlignedVector3(m_coeffs / norm()); } - + inline Scalar sum() const { ei_assert(m_coeffs.w()==Scalar(0)); return m_coeffs.sum(); } - + inline Scalar squaredNorm() const { ei_assert(m_coeffs.w()==Scalar(0)); return m_coeffs.squaredNorm(); } - + inline Scalar norm() const { return ei_sqrt(squaredNorm()); } - + inline AlignedVector3 cross(const AlignedVector3& other) const { return AlignedVector3(m_coeffs.cross3(other.m_coeffs)); } - + template<typename Derived> inline bool isApprox(const MatrixBase<Derived>& other, RealScalar eps=precision<Scalar>()) const {
diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h b/unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h index 0feb577..d421973 100644 --- a/unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h
@@ -56,7 +56,7 @@ typedef Matrix<ActiveScalar, InputsAtCompileTime, 1> ActiveInput; typedef Matrix<ActiveScalar, ValuesAtCompileTime, 1> ActiveValue; - void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const + void operator() (const InputType& x, ValueType* v, JacobianType* _jac=0) const { ei_assert(v!=0); if (!_jac)
diff --git a/unsupported/test/autodiff.cpp b/unsupported/test/autodiff.cpp index 43ef1a3..b116489 100644 --- a/unsupported/test/autodiff.cpp +++ b/unsupported/test/autodiff.cpp
@@ -111,7 +111,7 @@ } }; -template<typename Func> void adolc_forward_jacobian(const Func& f) +template<typename Func> void forward_jacobian(const Func& f) { typename Func::InputType x = Func::InputType::Random(f.inputs()); typename Func::ValueType y(f.values()), yref(f.values()); @@ -134,21 +134,29 @@ VERIFY_IS_APPROX(j, jref); } -void test_autodiff() +void test_autodiff_scalar() { std::cerr << foo<float>(1,2) << "\n"; AutoDiffScalar<Vector2f> ax(1,Vector2f::UnitX()); AutoDiffScalar<Vector2f> ay(2,Vector2f::UnitY()); std::cerr << foo<AutoDiffScalar<Vector2f> >(ax,ay).value() << " <> " << foo<AutoDiffScalar<Vector2f> >(ax,ay).derivatives().transpose() << "\n\n"; - - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,2>()) )); - CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,3>()) )); - CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,2>()) )); - CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) )); - CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double>(3,3)) )); - } - -// exit(1); } + +void test_autodiff_jacobian() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,2>()) )); + CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,3>()) )); + CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,2>()) )); + CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,3>()) )); + CALL_SUBTEST(( forward_jacobian(TestFunc1<double>(3,3)) )); + } +} + +void test_autodiff() +{ + test_autodiff_scalar(); + test_autodiff_jacobian(); +} +