From ac1bb3e5b3349bf744c3905fec55295e6bb024ab Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 18 Jul 2014 14:22:33 +0200 Subject: bug #770: fix out of bounds access --- unsupported/Eigen/MPRealSupport | 3 +++ 1 file changed, 3 insertions(+) diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport index 35d77e5bd..8e22b13c7 100644 --- a/unsupported/Eigen/MPRealSupport +++ b/unsupported/Eigen/MPRealSupport @@ -157,6 +157,9 @@ int main() void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha, Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0) { + if(rows==0 || cols==0 || depth==0) + return; + mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())), tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr())); -- cgit v1.2.3 From 1cb71a8782c5aa8d75562ebf2c36d1581e7799ef Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Fri, 18 Jul 2014 14:34:58 +0200 Subject: bug #138: Make building of internal documentation configurable via cmake flag --- doc/CMakeLists.txt | 8 ++++++++ doc/Doxyfile.in | 10 +++++----- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 8a493031c..1b8aaf9aa 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -10,12 +10,20 @@ if(CMAKE_COMPILER_IS_GNUCXX) endif(CMAKE_SYSTEM_NAME MATCHES Linux) endif(CMAKE_COMPILER_IS_GNUCXX) +option(EIGEN_INTERNAL_DOCUMENTATION "Build internal documentation" OFF) + + # Set some Doxygen flags set(EIGEN_DOXY_PROJECT_NAME "Eigen") set(EIGEN_DOXY_OUTPUT_DIRECTORY_SUFFIX "") set(EIGEN_DOXY_INPUT "\"${Eigen_SOURCE_DIR}/Eigen\" \"${Eigen_SOURCE_DIR}/doc\"") set(EIGEN_DOXY_HTML_COLORSTYLE_HUE "220") set(EIGEN_DOXY_TAGFILES "") +if(EIGEN_INTERNAL_DOCUMENTATION) + set(EIGEN_DOXY_INTERNAL "YES") +else(EIGEN_INTERNAL_DOCUMENTATION) + set(EIGEN_DOXY_INTERNAL "NO") +endif(EIGEN_INTERNAL_DOCUMENTATION) configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index 2ceee5e20..5d82add72 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -460,7 +460,7 @@ HIDE_IN_BODY_DOCS = NO # to NO (the default) then the documentation will be excluded. # Set it to YES to include the internal documentation. -INTERNAL_DOCS = NO +INTERNAL_DOCS = ${EIGEN_DOXY_INTERNAL} # If the CASE_SENSE_NAMES tag is set to NO then Doxygen will only generate # file names in lower-case letters. If set to YES upper-case letters are also @@ -480,7 +480,7 @@ HIDE_SCOPE_NAMES = NO # will put a list of the files that are included by a file in the documentation # of that file. -SHOW_INCLUDE_FILES = NO +SHOW_INCLUDE_FILES = ${EIGEN_DOXY_INTERNAL} # If the FORCE_LOCAL_INCLUDES tag is set to YES then Doxygen # will list include files with double quotes in the documentation @@ -546,7 +546,7 @@ STRICT_PROTO_MATCHING = NO # disable (NO) the todo list. This list is created by putting \todo # commands in the documentation. -GENERATE_TODOLIST = NO +GENERATE_TODOLIST = ${EIGEN_DOXY_INTERNAL} # The GENERATE_TESTLIST tag can be used to enable (YES) or # disable (NO) the test list. This list is created by putting \test @@ -558,13 +558,13 @@ GENERATE_TESTLIST = NO # disable (NO) the bug list. This list is created by putting \bug # commands in the documentation. -GENERATE_BUGLIST = NO +GENERATE_BUGLIST = ${EIGEN_DOXY_INTERNAL} # The GENERATE_DEPRECATEDLIST tag can be used to enable (YES) or # disable (NO) the deprecated list. This list is created by putting # \deprecated commands in the documentation. -GENERATE_DEPRECATEDLIST= NO +GENERATE_DEPRECATEDLIST= ${EIGEN_DOXY_INTERNAL} # The ENABLED_SECTIONS tag can be used to enable conditional # documentation sections, marked by \if sectionname ... \endif. -- cgit v1.2.3 From 68eafc10b1165eb83461187e549a54a618495916 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Fri, 18 Jul 2014 15:42:12 +0200 Subject: Add note to EIGEN_DONT_PARALLELIZE into preprocessor documentation page (requested in IRC) --- doc/PreprocessorDirectives.dox | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/PreprocessorDirectives.dox b/doc/PreprocessorDirectives.dox index e67991fb7..96d273823 100644 --- a/doc/PreprocessorDirectives.dox +++ b/doc/PreprocessorDirectives.dox @@ -61,6 +61,8 @@ run time. However, these assertions do cost time and can thus be turned off. expect that any objects passed to it are aligned. This will turn off vectorization. Not defined by default. - \b EIGEN_DONT_ALIGN_STATICALLY - disables alignment of arrays on the stack. Not defined by default, unless \c EIGEN_DONT_ALIGN is defined. + - \b EIGEN_DONT_PARALLELIZE - if defined, this disables multi-threading. This is only relevant if you enabled OpenMP. + See \ref TopicMultiThreading for details. - \b EIGEN_DONT_VECTORIZE - disables explicit vectorization when defined. Not defined by default, unless alignment is disabled by %Eigen's platform test or the user defining \c EIGEN_DONT_ALIGN. - \b EIGEN_FAST_MATH - enables some optimizations which might affect the accuracy of the result. This currently -- cgit v1.2.3 From ef4a86d6b8eb7334ffa3a5a55f8432d723f8a502 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Fri, 18 Jul 2014 16:39:58 +0200 Subject: Fix trivial warnings in MPRealSupport --- unsupported/Eigen/MPRealSupport | 6 +++--- unsupported/test/mpreal/mpreal.h | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport index 8e22b13c7..632de3854 100644 --- a/unsupported/Eigen/MPRealSupport +++ b/unsupported/Eigen/MPRealSupport @@ -88,9 +88,9 @@ int main() inline static Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); } inline static Real dummy_precision() - { - unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; - return mpfr::machine_epsilon(weak_prec); + { + mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; + return mpfr::machine_epsilon(weak_prec); } }; diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index 4c29d2ac2..dddda7dd9 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -1698,7 +1698,7 @@ inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcpt ////////////////////////////////////////////////////////////////////////// // Type Converters -inline bool mpreal::toBool (mp_rnd_t mode) const { return mpfr_zero_p (mpfr_srcptr()) == 0; } +inline bool mpreal::toBool (mp_rnd_t /*mode*/) const { return mpfr_zero_p (mpfr_srcptr()) == 0; } inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si (mpfr_srcptr(), mode); } inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui (mpfr_srcptr(), mode); } inline float mpreal::toFloat (mp_rnd_t mode) const { return mpfr_get_flt(mpfr_srcptr(), mode); } @@ -3070,4 +3070,4 @@ namespace std } -#endif /* __MPREAL_H__ */ \ No newline at end of file +#endif /* __MPREAL_H__ */ -- cgit v1.2.3 From 8e19027130b1e31cfab4e09b67c155c07cf85ff6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 20 Jul 2014 14:03:22 +0200 Subject: bug #826: fix 64 to 32 bits conversion warning when calling Matrix(long) --- Eigen/src/Core/PlainObjectBase.h | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/Eigen/src/Core/PlainObjectBase.h b/Eigen/src/Core/PlainObjectBase.h index ae5b342ce..e4e3e8903 100644 --- a/Eigen/src/Core/PlainObjectBase.h +++ b/Eigen/src/Core/PlainObjectBase.h @@ -705,6 +705,17 @@ class PlainObjectBase : public internal::dense_xpr_base::type EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 1) m_storage.data()[0] = val0; } + + template + EIGEN_DEVICE_FUNC + EIGEN_STRONG_INLINE void _init1(const Index& val0, + typename internal::enable_if< (!internal::is_same::value) + && Base::SizeAtCompileTime==1 + && internal::is_convertible::value,T>::type* = 0) + { + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 1) + m_storage.data()[0] = Scalar(val0); + } template EIGEN_DEVICE_FUNC -- cgit v1.2.3 From d4cc1bdc7f27b856c4aefc939254ec851c831bf0 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 20 Jul 2014 14:22:58 +0200 Subject: Make the ordering method of SimplicialL[D]LT user configurable. --- Eigen/src/SparseCholesky/SimplicialCholesky.h | 42 +++++++++++++++------------ test/simplicial_cholesky.cpp | 41 ++++++++++++++------------ 2 files changed, 46 insertions(+), 37 deletions(-) diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index f41d7e010..e1f96ba5a 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -37,6 +37,7 @@ class SimplicialCholeskyBase : internal::noncopyable { public: typedef typename internal::traits::MatrixType MatrixType; + typedef typename internal::traits::OrderingType OrderingType; enum { UpLo = internal::traits::UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -240,15 +241,16 @@ class SimplicialCholeskyBase : internal::noncopyable RealScalar m_shiftScale; }; -template class SimplicialLLT; -template class SimplicialLDLT; -template class SimplicialCholesky; +template > class SimplicialLLT; +template > class SimplicialLDLT; +template > class SimplicialCholesky; namespace internal { -template struct traits > +template struct traits > { typedef _MatrixType MatrixType; + typedef _Ordering OrderingType; enum { UpLo = _UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -259,9 +261,10 @@ template struct traits struct traits > +template struct traits > { typedef _MatrixType MatrixType; + typedef _Ordering OrderingType; enum { UpLo = _UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -272,9 +275,10 @@ template struct traits struct traits > +template struct traits > { typedef _MatrixType MatrixType; + typedef _Ordering OrderingType; enum { UpLo = _UpLo }; }; @@ -294,11 +298,12 @@ template struct traits * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower * or Upper. Default is Lower. + * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> * - * \sa class SimplicialLDLT + * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering */ -template - class SimplicialLLT : public SimplicialCholeskyBase > +template + class SimplicialLLT : public SimplicialCholeskyBase > { public: typedef _MatrixType MatrixType; @@ -382,11 +387,12 @@ public: * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower * or Upper. Default is Lower. + * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> * - * \sa class SimplicialLLT + * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering */ -template - class SimplicialLDLT : public SimplicialCholeskyBase > +template + class SimplicialLDLT : public SimplicialCholeskyBase > { public: typedef _MatrixType MatrixType; @@ -467,8 +473,8 @@ public: * * \sa class SimplicialLDLT, class SimplicialLLT */ -template - class SimplicialCholesky : public SimplicialCholeskyBase > +template + class SimplicialCholesky : public SimplicialCholeskyBase > { public: typedef _MatrixType MatrixType; @@ -612,15 +618,13 @@ void SimplicialCholeskyBase::ordering(const MatrixType& a, CholMatrixTy { eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); - // TODO allows to configure the permutation // Note that amd compute the inverse permutation { CholMatrixType C; C = a.template selfadjointView(); - // remove diagonal entries: - // seems not to be needed - // C.prune(keep_diag()); - internal::minimum_degree_ordering(C, m_Pinv); + + OrderingType ordering; + ordering(C,m_Pinv); } if(m_Pinv.size()>0) diff --git a/test/simplicial_cholesky.cpp b/test/simplicial_cholesky.cpp index e93a52e9c..786468421 100644 --- a/test/simplicial_cholesky.cpp +++ b/test/simplicial_cholesky.cpp @@ -11,26 +11,31 @@ template void test_simplicial_cholesky_T() { - SimplicialCholesky, Lower> chol_colmajor_lower; - SimplicialCholesky, Upper> chol_colmajor_upper; - SimplicialLLT, Lower> llt_colmajor_lower; - SimplicialLDLT, Upper> llt_colmajor_upper; - SimplicialLDLT, Lower> ldlt_colmajor_lower; - SimplicialLDLT, Upper> ldlt_colmajor_upper; + SimplicialCholesky, Lower> chol_colmajor_lower_amd; + SimplicialCholesky, Upper> chol_colmajor_upper_amd; + SimplicialLLT, Lower> llt_colmajor_lower_amd; + SimplicialLLT, Upper> llt_colmajor_upper_amd; + SimplicialLDLT, Lower> ldlt_colmajor_lower_amd; + SimplicialLDLT, Upper> ldlt_colmajor_upper_amd; + SimplicialLDLT, Lower, NaturalOrdering > ldlt_colmajor_lower_nat; + SimplicialLDLT, Upper, NaturalOrdering > ldlt_colmajor_upper_nat; - check_sparse_spd_solving(chol_colmajor_lower); - check_sparse_spd_solving(chol_colmajor_upper); - check_sparse_spd_solving(llt_colmajor_lower); - check_sparse_spd_solving(llt_colmajor_upper); - check_sparse_spd_solving(ldlt_colmajor_lower); - check_sparse_spd_solving(ldlt_colmajor_upper); + check_sparse_spd_solving(chol_colmajor_lower_amd); + check_sparse_spd_solving(chol_colmajor_upper_amd); + check_sparse_spd_solving(llt_colmajor_lower_amd); + check_sparse_spd_solving(llt_colmajor_upper_amd); + check_sparse_spd_solving(ldlt_colmajor_lower_amd); + check_sparse_spd_solving(ldlt_colmajor_upper_amd); - check_sparse_spd_determinant(chol_colmajor_lower); - check_sparse_spd_determinant(chol_colmajor_upper); - check_sparse_spd_determinant(llt_colmajor_lower); - check_sparse_spd_determinant(llt_colmajor_upper); - check_sparse_spd_determinant(ldlt_colmajor_lower); - check_sparse_spd_determinant(ldlt_colmajor_upper); + check_sparse_spd_determinant(chol_colmajor_lower_amd); + check_sparse_spd_determinant(chol_colmajor_upper_amd); + check_sparse_spd_determinant(llt_colmajor_lower_amd); + check_sparse_spd_determinant(llt_colmajor_upper_amd); + check_sparse_spd_determinant(ldlt_colmajor_lower_amd); + check_sparse_spd_determinant(ldlt_colmajor_upper_amd); + + check_sparse_spd_solving(ldlt_colmajor_lower_nat); + check_sparse_spd_solving(ldlt_colmajor_upper_nat); } void test_simplicial_cholesky() -- cgit v1.2.3 From 339f14b8d1b73db0366afc3a497d566cecf72e1b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 21 Jul 2014 13:43:48 +0200 Subject: bug #826: document caveats in 1x1 and 2x1 constructors. --- Eigen/src/Core/Matrix.h | 22 ++++++++++++++++++---- doc/PreprocessorDirectives.dox | 13 +++++++++++++ 2 files changed, 31 insertions(+), 4 deletions(-) diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index cd2ff91e7..8c95ee3ca 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -257,9 +257,15 @@ class Matrix /** \brief Constructs a vector or row-vector with given dimension. \only_for_vectors * - * Note that this is only useful for dynamic-size vectors. For fixed-size vectors, - * it is redundant to pass the dimension here, so it makes more sense to use the default - * constructor Matrix() instead. + * This is useful for dynamic-size vectors. For fixed-size vectors, + * it is redundant to pass these parameters, so one should use the default constructor + * Matrix() instead. + * + * \warning This constructor is disabled for fixed-size \c 1x1 matrices. For instance, + * calling Matrix(1) will call the initialization constructor: Matrix(const Scalar&). + * For fixed-size \c 1x1 matrices it is thefore recommended to use the default + * constructor Matrix() instead, especilly when using one of the non standard + * \c EIGEN_INITIALIZE_MATRICES_BY_{ZERO,\c NAN} macros (see \ref TopicPreprocessorDirectives). */ EIGEN_STRONG_INLINE explicit Matrix(Index dim); /** \brief Constructs an initialized 1x1 matrix with the given coefficient */ @@ -268,9 +274,17 @@ class Matrix * * This is useful for dynamic-size matrices. For fixed-size matrices, * it is redundant to pass these parameters, so one should use the default constructor - * Matrix() instead. */ + * Matrix() instead. + * + * \warning This constructor is disabled for fixed-size \c 1x2 and \c 2x1 vectors. For instance, + * calling Matrix2f(2,1) will call the initialization constructor: Matrix(const Scalar& x, const Scalar& y). + * For fixed-size \c 1x2 or \c 2x1 vectors it is thefore recommended to use the default + * constructor Matrix() instead, especilly when using one of the non standard + * \c EIGEN_INITIALIZE_MATRICES_BY_{ZERO,\c NAN} macros (see \ref TopicPreprocessorDirectives). + */ EIGEN_DEVICE_FUNC Matrix(Index rows, Index cols); + /** \brief Constructs an initialized 2D vector with given coefficients */ Matrix(const Scalar& x, const Scalar& y); #endif diff --git a/doc/PreprocessorDirectives.dox b/doc/PreprocessorDirectives.dox index 96d273823..4be2167ef 100644 --- a/doc/PreprocessorDirectives.dox +++ b/doc/PreprocessorDirectives.dox @@ -27,10 +27,23 @@ are doing. Defaults to the %IOFormat constructed by the default constructor IOFormat::IOFormat(). - \b EIGEN_INITIALIZE_MATRICES_BY_ZERO - if defined, all entries of newly constructed matrices and arrays are initialized to zero, as are new entries in matrices and arrays after resizing. Not defined by default. + \warning The unary (resp. binary) constructor of \c 1x1 (resp. \c 2x1 or \c 1x2) fixed size matrices is + always interpreted as an initialization constructor where the argument(s) are the coefficient values + and not the sizes. For instance, \code Vector2d v(2,1); \endcode will create a vector with coeficients [2,1], + and \b not a \c 2x1 vector initialized with zeros (i.e., [0,0]). If such cases might occur, then it is + recommended to use the default constructor with a explicit call to resize: + \code + Matrix v; + v.resize(size); + Matrix m; + m.resize(rows,cols); + \endcode - \b EIGEN_INITIALIZE_MATRICES_BY_NAN - if defined, all entries of newly constructed matrices and arrays are initialized to NaN, as are new entries in matrices and arrays after resizing. This option is especially useful for debugging purpose, though a memory tool like valgrind is preferable. Not defined by default. + \warning See the documentation of \c EIGEN_INITIALIZE_MATRICES_BY_ZERO for a discussion on a limitations + of these macros when applied to \c 1x1, \c 1x2, and \c 2x1 fixed-size matrices. - \b EIGEN_NO_AUTOMATIC_RESIZING - if defined, the matrices (or arrays) on both sides of an assignment a = b have to be of the same size; otherwise, %Eigen automatically resizes \c a so that it is of the correct size. Not defined by default. -- cgit v1.2.3 From 58687aa5e638d365d5e41c1e6c66cbfc44fce85f Mon Sep 17 00:00:00 2001 From: Moritz Klammler Date: Sun, 6 Jul 2014 06:58:13 +0200 Subject: Avoid memory leak when constructor of user-defined type throws exception. The added check `ctorleak.cpp` demonstrates how the leak can be reproduced. The test appears to pass but it is leaking the storage of the (not created) matrix. I don't know how to make this test fail in the existing test suite but you can run it through Valgrind (or another debugger) to verify the leak. $ ./check.sh ctorleak && valgrind --leak-check=full ./test/ctorleak This patch fixes this leak by adding some try-catch-delete-rethrow blocks to `Eigen/src/Core/util/Memory.h`. --- Eigen/src/Core/util/Memory.h | 83 ++++++++++++++++++++++++++++++++++++-------- test/CMakeLists.txt | 2 ++ test/ctorleak.cpp | 43 +++++++++++++++++++++++ 3 files changed, 114 insertions(+), 14 deletions(-) create mode 100644 test/ctorleak.cpp diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index dd9285714..4b3dcde3a 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -338,15 +338,6 @@ template<> inline void* conditional_aligned_realloc(void* ptr, size_t new *** Construction/destruction of array elements *** *****************************************************************************/ -/** \internal Constructs the elements of an array. - * The \a size parameter tells on how many objects to call the constructor of T. - */ -template inline T* construct_elements_of_array(T *ptr, size_t size) -{ - for (size_t i=0; i < size; ++i) ::new (ptr + i) T; - return ptr; -} - /** \internal Destructs the elements of an array. * The \a size parameters tells on how many objects to call the destructor of T. */ @@ -357,6 +348,24 @@ template inline void destruct_elements_of_array(T *ptr, size_t size) while(size) ptr[--size].~T(); } +/** \internal Constructs the elements of an array. + * The \a size parameter tells on how many objects to call the constructor of T. + */ +template inline T* construct_elements_of_array(T *ptr, size_t size) +{ + size_t i; + try + { + for (i = 0; i < size; ++i) ::new (ptr + i) T; + return ptr; + } + catch (...) + { + destruct_elements_of_array(ptr, i); + throw; + } +} + /***************************************************************************** *** Implementation of aligned new/delete-like functions *** *****************************************************************************/ @@ -376,14 +385,30 @@ template inline T* aligned_new(size_t size) { check_size_for_overflow(size); T *result = reinterpret_cast(aligned_malloc(sizeof(T)*size)); - return construct_elements_of_array(result, size); + try + { + return construct_elements_of_array(result, size); + } + catch (...) + { + aligned_free(result); + throw; + } } template inline T* conditional_aligned_new(size_t size) { check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); - return construct_elements_of_array(result, size); + try + { + return construct_elements_of_array(result, size); + } + catch (...) + { + conditional_aligned_free(result); + throw; + } } /** \internal Deletes objects constructed with aligned_new @@ -412,7 +437,17 @@ template inline T* conditional_aligned_realloc_new(T* pt destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); if(new_size > old_size) - construct_elements_of_array(result+old_size, new_size-old_size); + { + try + { + construct_elements_of_array(result+old_size, new_size-old_size); + } + catch (...) + { + conditional_aligned_free(result); + throw; + } + } return result; } @@ -422,7 +457,17 @@ template inline T* conditional_aligned_new_auto(size_t s check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); if(NumTraits::RequireInitialization) - construct_elements_of_array(result, size); + { + try + { + construct_elements_of_array(result, size); + } + catch (...) + { + conditional_aligned_free(result); + throw; + } + } return result; } @@ -434,7 +479,17 @@ template inline T* conditional_aligned_realloc_new_auto( destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); if(NumTraits::RequireInitialization && (new_size > old_size)) - construct_elements_of_array(result+old_size, new_size-old_size); + { + try + { + construct_elements_of_array(result+old_size, new_size-old_size); + } + catch (...) + { + conditional_aligned_free(result); + throw; + } + } return result; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4521d07e4..fed5c0e06 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -251,6 +251,8 @@ ei_add_test(bicgstab) ei_add_test(sparselu) ei_add_test(sparseqr) +ei_add_test(ctorleak) + # ei_add_test(denseLM) if(QT4_FOUND) diff --git a/test/ctorleak.cpp b/test/ctorleak.cpp new file mode 100644 index 000000000..72ab94b66 --- /dev/null +++ b/test/ctorleak.cpp @@ -0,0 +1,43 @@ +#include "main.h" + +#include // std::exception + +struct Foo +{ + int dummy; + Foo() { if (!internal::random(0, 10)) throw Foo::Fail(); } + class Fail : public std::exception {}; +}; + +namespace Eigen +{ + template<> + struct NumTraits + { + typedef double Real; + typedef double NonInteger; + typedef double Nested; + enum + { + IsComplex = 0, + IsInteger = 1, + ReadCost = -1, + AddCost = -1, + MulCost = -1, + IsSigned = 1, + RequireInitialization = 1 + }; + static inline Real epsilon() { return 1.0; } + static inline Real dummy_epsilon() { return 0.0; } + }; +} + +void test_ctorleak() +{ + try + { + Matrix m(14, 92); + eigen_assert(false); // not reached + } + catch (const Foo::Fail&) { /* ignore */ } +} -- cgit v1.2.3 From 529e6cb5529f626b7c5cf781bdcd5a5720f904c1 Mon Sep 17 00:00:00 2001 From: Moritz Klammler Date: Fri, 18 Jul 2014 23:19:56 +0200 Subject: Applied changes suggested by Christoph Hertzberg to c'tor leak fix. - Enclose exception handling in '#ifdef EIGEN_EXCEPTIONS'. - Use an object counter to demonstrate the bug more readily. --- Eigen/src/Core/util/Memory.h | 24 ++++++++++++++++++++++++ test/ctorleak.cpp | 28 +++++++++++++++++++++++++++- 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 4b3dcde3a..9b3e84fb0 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -354,16 +354,20 @@ template inline void destruct_elements_of_array(T *ptr, size_t size) template inline T* construct_elements_of_array(T *ptr, size_t size) { size_t i; +#ifdef EIGEN_EXCEPTIONS try +#endif { for (i = 0; i < size; ++i) ::new (ptr + i) T; return ptr; } +#ifdef EIGEN_EXCEPTIONS catch (...) { destruct_elements_of_array(ptr, i); throw; } +#endif } /***************************************************************************** @@ -385,30 +389,38 @@ template inline T* aligned_new(size_t size) { check_size_for_overflow(size); T *result = reinterpret_cast(aligned_malloc(sizeof(T)*size)); +#ifdef EIGEN_EXCEPTIONS try +#endif { return construct_elements_of_array(result, size); } +#ifdef EIGEN_EXCEPTIONS catch (...) { aligned_free(result); throw; } +#endif } template inline T* conditional_aligned_new(size_t size) { check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); +#ifdef EIGEN_EXCEPTIONS try +#endif { return construct_elements_of_array(result, size); } +#ifdef EIGEN_EXCEPTIONS catch (...) { conditional_aligned_free(result); throw; } +#endif } /** \internal Deletes objects constructed with aligned_new @@ -438,15 +450,19 @@ template inline T* conditional_aligned_realloc_new(T* pt T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); if(new_size > old_size) { +#ifdef EIGEN_EXCEPTIONS try +#endif { construct_elements_of_array(result+old_size, new_size-old_size); } +#ifdef EIGEN_EXCEPTIONS catch (...) { conditional_aligned_free(result); throw; } +#endif } return result; } @@ -458,15 +474,19 @@ template inline T* conditional_aligned_new_auto(size_t s T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); if(NumTraits::RequireInitialization) { +#ifdef EIGEN_EXCEPTIONS try +#endif { construct_elements_of_array(result, size); } +#ifdef EIGEN_EXCEPTIONS catch (...) { conditional_aligned_free(result); throw; } +#endif } return result; } @@ -480,15 +500,19 @@ template inline T* conditional_aligned_realloc_new_auto( T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); if(NumTraits::RequireInitialization && (new_size > old_size)) { +#ifdef EIGEN_EXCEPTIONS try +#endif { construct_elements_of_array(result+old_size, new_size-old_size); } +#ifdef EIGEN_EXCEPTIONS catch (...) { conditional_aligned_free(result); throw; } +#endif } return result; } diff --git a/test/ctorleak.cpp b/test/ctorleak.cpp index 72ab94b66..f3f4411c8 100644 --- a/test/ctorleak.cpp +++ b/test/ctorleak.cpp @@ -4,11 +4,30 @@ struct Foo { + static unsigned object_count; + static unsigned object_limit; int dummy; - Foo() { if (!internal::random(0, 10)) throw Foo::Fail(); } + + Foo() + { +#ifdef EIGEN_EXCEPTIONS + // TODO: Is this the correct way to handle this? + if (Foo::object_count > Foo::object_limit) { throw Foo::Fail(); } +#endif + ++Foo::object_count; + } + + ~Foo() + { + --Foo::object_count; + } + class Fail : public std::exception {}; }; +unsigned Foo::object_count = 0; +unsigned Foo::object_limit = 0; + namespace Eigen { template<> @@ -34,10 +53,17 @@ namespace Eigen void test_ctorleak() { + Foo::object_count = 0; + Foo::object_limit = internal::random(0, 14 * 92 - 2); +#ifdef EIGEN_EXCEPTIONS try +#endif { Matrix m(14, 92); eigen_assert(false); // not reached } +#ifdef EIGEN_EXCEPTIONS catch (const Foo::Fail&) { /* ignore */ } +#endif + VERIFY_IS_EQUAL(static_cast(0), Foo::object_count); } -- cgit v1.2.3 From a8283e0ed2d3e42ba8e1c42f51cb08eacc301047 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Tue, 22 Jul 2014 13:16:44 +0200 Subject: Define EIGEN_TRY, EIGEN_CATCH, EIGEN_THROW as suggested by Moritz Klammer. Make it possible to run unit-tests with exceptions disabled via EIGEN_TEST_NO_EXCEPTIONS flag. Enhanced ctorleak unit-test --- CMakeLists.txt | 8 ++- Eigen/Core | 16 ++--- Eigen/src/Core/util/Macros.h | 12 ++++ Eigen/src/Core/util/Memory.h | 145 +++++++++++++++++-------------------------- test/CMakeLists.txt | 10 +-- test/ctorleak.cpp | 40 ++++-------- test/main.h | 16 +++-- 7 files changed, 112 insertions(+), 135 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 470095680..96d6c8701 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -119,7 +119,7 @@ endmacro(ei_add_cxx_compiler_flag) if(NOT MSVC) # We assume that other compilers are partly compatible with GNUCC - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fexceptions") +# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fexceptions") set(CMAKE_CXX_FLAGS_DEBUG "-g3") set(CMAKE_CXX_FLAGS_RELEASE "-g0 -O2") @@ -301,6 +301,12 @@ if(EIGEN_TEST_NO_EXPLICIT_ALIGNMENT) message(STATUS "Disabling alignment in tests/examples") endif() +option(EIGEN_TEST_NO_EXCEPTIONS "Disables C++ exceptions" OFF) +if(EIGEN_TEST_NO_EXCEPTIONS) + ei_add_cxx_compiler_flag("-fno-exceptions") + message(STATUS "Disabling exceptions in tests/examples") +endif() + option(EIGEN_TEST_C++0x "Enables all C++0x features." OFF) include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/Eigen/Core b/Eigen/Core index 16e388fa4..9a73fe37b 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -42,6 +42,14 @@ #define EIGEN_USING_STD_MATH(FUNC) using std::FUNC; #endif +#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(__CUDA_ARCH__) + #define EIGEN_EXCEPTIONS +#endif + +#ifdef EIGEN_EXCEPTIONS + #include +#endif + // then include this file where all our macros are defined. It's really important to do it first because // it's where we do all the alignment settings (platform detection and honoring the user's will if he // defined e.g. EIGEN_DONT_ALIGN) so it needs to be done before we do anything with vectorization. @@ -209,14 +217,6 @@ #include #endif -#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(__CUDA_ARCH__) - #define EIGEN_EXCEPTIONS -#endif - -#ifdef EIGEN_EXCEPTIONS - #include -#endif - /** \brief Namespace containing all symbols from the %Eigen library. */ namespace Eigen { diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 71e42b125..5e9b0a112 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -452,4 +452,16 @@ namespace Eigen { const RHS \ > +#ifdef EIGEN_EXCEPTIONS +# define EIGEN_THROW_X(X) throw X +# define EIGEN_THROW throw +# define EIGEN_TRY try +# define EIGEN_CATCH(X) catch (X) +#else +# define EIGEN_THROW_X(X) std::abort() +# define EIGEN_THROW std::abort() +# define EIGEN_TRY if (true) +# define EIGEN_CATCH(X) else +#endif + #endif // EIGEN_MACROS_H diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 9b3e84fb0..38990566d 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -354,20 +354,16 @@ template inline void destruct_elements_of_array(T *ptr, size_t size) template inline T* construct_elements_of_array(T *ptr, size_t size) { size_t i; -#ifdef EIGEN_EXCEPTIONS - try -#endif - { + EIGEN_TRY + { for (i = 0; i < size; ++i) ::new (ptr + i) T; return ptr; - } -#ifdef EIGEN_EXCEPTIONS - catch (...) - { - destruct_elements_of_array(ptr, i); - throw; - } -#endif + } + EIGEN_CATCH(...) + { + destruct_elements_of_array(ptr, i); + EIGEN_THROW; + } } /***************************************************************************** @@ -389,38 +385,30 @@ template inline T* aligned_new(size_t size) { check_size_for_overflow(size); T *result = reinterpret_cast(aligned_malloc(sizeof(T)*size)); -#ifdef EIGEN_EXCEPTIONS - try -#endif - { - return construct_elements_of_array(result, size); - } -#ifdef EIGEN_EXCEPTIONS - catch (...) - { - aligned_free(result); - throw; - } -#endif + EIGEN_TRY + { + return construct_elements_of_array(result, size); + } + EIGEN_CATCH(...) + { + aligned_free(result); + EIGEN_THROW; + } } template inline T* conditional_aligned_new(size_t size) { check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); -#ifdef EIGEN_EXCEPTIONS - try -#endif - { - return construct_elements_of_array(result, size); - } -#ifdef EIGEN_EXCEPTIONS - catch (...) - { - conditional_aligned_free(result); - throw; - } -#endif + EIGEN_TRY + { + return construct_elements_of_array(result, size); + } + EIGEN_CATCH(...) + { + conditional_aligned_free(result); + EIGEN_THROW; + } } /** \internal Deletes objects constructed with aligned_new @@ -449,21 +437,17 @@ template inline T* conditional_aligned_realloc_new(T* pt destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); if(new_size > old_size) + { + EIGEN_TRY { -#ifdef EIGEN_EXCEPTIONS - try -#endif - { - construct_elements_of_array(result+old_size, new_size-old_size); - } -#ifdef EIGEN_EXCEPTIONS - catch (...) - { - conditional_aligned_free(result); - throw; - } -#endif + construct_elements_of_array(result+old_size, new_size-old_size); + } + EIGEN_CATCH(...) + { + conditional_aligned_free(result); + EIGEN_THROW; } + } return result; } @@ -473,21 +457,17 @@ template inline T* conditional_aligned_new_auto(size_t s check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); if(NumTraits::RequireInitialization) + { + EIGEN_TRY { -#ifdef EIGEN_EXCEPTIONS - try -#endif - { - construct_elements_of_array(result, size); - } -#ifdef EIGEN_EXCEPTIONS - catch (...) - { - conditional_aligned_free(result); - throw; - } -#endif + construct_elements_of_array(result, size); } + EIGEN_CATCH(...) + { + conditional_aligned_free(result); + EIGEN_THROW; + } + } return result; } @@ -499,21 +479,17 @@ template inline T* conditional_aligned_realloc_new_auto( destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); if(NumTraits::RequireInitialization && (new_size > old_size)) + { + EIGEN_TRY { -#ifdef EIGEN_EXCEPTIONS - try -#endif - { - construct_elements_of_array(result+old_size, new_size-old_size); - } -#ifdef EIGEN_EXCEPTIONS - catch (...) - { - conditional_aligned_free(result); - throw; - } -#endif + construct_elements_of_array(result+old_size, new_size-old_size); + } + EIGEN_CATCH(...) + { + conditional_aligned_free(result); + EIGEN_THROW; } + } return result; } @@ -713,20 +689,11 @@ template class aligned_stack_memory_handler *****************************************************************************/ #if EIGEN_ALIGN - #ifdef EIGEN_EXCEPTIONS - #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \ - void* operator new(size_t size, const std::nothrow_t&) throw() { \ - try { return Eigen::internal::conditional_aligned_malloc(size); } \ - catch (...) { return 0; } \ - return 0; \ - } - #else - #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \ + #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \ void* operator new(size_t size, const std::nothrow_t&) throw() { \ - return Eigen::internal::conditional_aligned_malloc(size); \ + EIGEN_TRY { return Eigen::internal::conditional_aligned_malloc(size); } \ + EIGEN_CATCH (...) { return 0; } \ } - #endif - #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) \ void *operator new(size_t size) { \ return Eigen::internal::conditional_aligned_malloc(size); \ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fed5c0e06..47aefddb8 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -158,7 +158,9 @@ ei_add_test(basicstuff) ei_add_test(linearstructure) ei_add_test(integer_types) ei_add_test(unalignedcount) -ei_add_test(exceptions) +if(NOT EIGEN_TEST_NO_EXCEPTIONS) + ei_add_test(exceptions) +endif() ei_add_test(redux) ei_add_test(visitor) ei_add_test(block) @@ -238,7 +240,9 @@ ei_add_test(nesting_ops "${CMAKE_CXX_FLAGS_DEBUG}") ei_add_test(zerosized) ei_add_test(dontalign) ei_add_test(evaluators) -ei_add_test(sizeoverflow) +if(NOT EIGEN_TEST_NO_EXCEPTIONS) + ei_add_test(sizeoverflow) +endif() ei_add_test(prec_inverse_4x4) ei_add_test(vectorwiseop) ei_add_test(special_numbers) @@ -251,8 +255,6 @@ ei_add_test(bicgstab) ei_add_test(sparselu) ei_add_test(sparseqr) -ei_add_test(ctorleak) - # ei_add_test(denseLM) if(QT4_FOUND) diff --git a/test/ctorleak.cpp b/test/ctorleak.cpp index f3f4411c8..145d91be4 100644 --- a/test/ctorleak.cpp +++ b/test/ctorleak.cpp @@ -28,42 +28,24 @@ struct Foo unsigned Foo::object_count = 0; unsigned Foo::object_limit = 0; -namespace Eigen -{ - template<> - struct NumTraits - { - typedef double Real; - typedef double NonInteger; - typedef double Nested; - enum - { - IsComplex = 0, - IsInteger = 1, - ReadCost = -1, - AddCost = -1, - MulCost = -1, - IsSigned = 1, - RequireInitialization = 1 - }; - static inline Real epsilon() { return 1.0; } - static inline Real dummy_epsilon() { return 0.0; } - }; -} void test_ctorleak() { + typedef DenseIndex Index; Foo::object_count = 0; - Foo::object_limit = internal::random(0, 14 * 92 - 2); + for(int i = 0; i < g_repeat; i++) { + Index rows = internal::random(2,EIGEN_TEST_MAX_SIZE), cols = internal::random(2,EIGEN_TEST_MAX_SIZE); + Foo::object_limit = internal::random(0, rows*cols - 2); #ifdef EIGEN_EXCEPTIONS - try -#endif + try { - Matrix m(14, 92); - eigen_assert(false); // not reached - } +#endif + Matrix m(rows, cols); #ifdef EIGEN_EXCEPTIONS - catch (const Foo::Fail&) { /* ignore */ } + VERIFY(false); // not reached if exceptions are enabled + } + catch (const Foo::Fail&) { /* ignore */ } #endif + } VERIFY_IS_EQUAL(static_cast(0), Foo::object_count); } diff --git a/test/main.h b/test/main.h index 3ccc2ae88..3295dcb71 100644 --- a/test/main.h +++ b/test/main.h @@ -117,13 +117,14 @@ namespace Eigen if(report_on_cerr_on_assert_failure) \ std::cerr << #a << " " __FILE__ << "(" << __LINE__ << ")\n"; \ Eigen::no_more_assert = true; \ - throw Eigen::eigen_assert_exception(); \ + EIGEN_THROW_X(Eigen::eigen_assert_exception()); \ } \ else if (Eigen::internal::push_assert) \ { \ eigen_assert_list.push_back(std::string(EI_PP_MAKE_STRING(__FILE__) " (" EI_PP_MAKE_STRING(__LINE__) ") : " #a) ); \ } + #ifdef EIGEN_EXCEPTIONS #define VERIFY_RAISES_ASSERT(a) \ { \ Eigen::no_more_assert = false; \ @@ -142,6 +143,7 @@ namespace Eigen Eigen::report_on_cerr_on_assert_failure = true; \ Eigen::internal::push_assert = false; \ } + #endif //EIGEN_EXCEPTIONS #elif !defined(__CUDACC__) // EIGEN_DEBUG_ASSERTS // see bug 89. The copy_bool here is working around a bug in gcc <= 4.3 @@ -152,9 +154,10 @@ namespace Eigen if(report_on_cerr_on_assert_failure) \ eigen_plain_assert(a); \ else \ - throw Eigen::eigen_assert_exception(); \ + EIGEN_THROW_X(Eigen::eigen_assert_exception()); \ } - #define VERIFY_RAISES_ASSERT(a) { \ + #ifdef EIGEN_EXCEPTIONS + #define VERIFY_RAISES_ASSERT(a) { \ Eigen::no_more_assert = false; \ Eigen::report_on_cerr_on_assert_failure = false; \ try { \ @@ -164,9 +167,14 @@ namespace Eigen catch (Eigen::eigen_assert_exception&) { VERIFY(true); } \ Eigen::report_on_cerr_on_assert_failure = true; \ } - + #endif //EIGEN_EXCEPTIONS #endif // EIGEN_DEBUG_ASSERTS +#ifndef VERIFY_RAISES_ASSERT + #define VERIFY_RAISES_ASSERT(a) \ + std::cout << "Can't VERIFY_RAISES_ASSERT( " #a " ) with exceptions disabled"; +#endif + #if !defined(__CUDACC__) #define EIGEN_USE_CUSTOM_ASSERT #endif -- cgit v1.2.3 From 922694a2d1531b6c50fc26a96b7d998089da72bf Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Jul 2014 16:57:14 +0200 Subject: Extend fixed-size ctor unit test and fix conversion warning. --- Eigen/src/Core/PlainObjectBase.h | 10 ++++++++++ test/basicstuff.cpp | 28 +++++++++++++++++++++++----- 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Core/PlainObjectBase.h b/Eigen/src/Core/PlainObjectBase.h index e4e3e8903..b968b325e 100644 --- a/Eigen/src/Core/PlainObjectBase.h +++ b/Eigen/src/Core/PlainObjectBase.h @@ -681,6 +681,7 @@ class PlainObjectBase : public internal::dense_xpr_base::type FLOATING_POINT_ARGUMENT_PASSED__INTEGER_WAS_EXPECTED) resize(nbRows,nbCols); } + template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void _init2(const Scalar& val0, const Scalar& val1, typename internal::enable_if::type* = 0) @@ -689,6 +690,15 @@ class PlainObjectBase : public internal::dense_xpr_base::type m_storage.data()[0] = val0; m_storage.data()[1] = val1; } + + template + EIGEN_DEVICE_FUNC + EIGEN_STRONG_INLINE void _init2(const Index& val0, const Index& val1, typename internal::enable_if<(!internal::is_same::value) && Base::SizeAtCompileTime==2,T1>::type* = 0) + { + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 2) + m_storage.data()[0] = Scalar(val0); + m_storage.data()[1] = Scalar(val1); + } template EIGEN_DEVICE_FUNC diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 56370d591..70cefe19f 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -183,6 +183,7 @@ void fixedSizeMatrixConstruction() Scalar raw[4]; for(int k=0; k<4; ++k) raw[k] = internal::random(); + { Matrix m(raw); Array a(raw); @@ -200,18 +201,34 @@ void fixedSizeMatrixConstruction() VERIFY((a==Array(raw[0],raw[1],raw[2])).all()); } { - Matrix m(raw); - Array a(raw); + Matrix m(raw), m2( (DenseIndex(raw[0])), (DenseIndex(raw[1])) ); + Array a(raw), a2( (DenseIndex(raw[0])), (DenseIndex(raw[1])) ); for(int k=0; k<2; ++k) VERIFY(m(k) == raw[k]); for(int k=0; k<2; ++k) VERIFY(a(k) == raw[k]); VERIFY_IS_EQUAL(m,(Matrix(raw[0],raw[1]))); VERIFY((a==Array(raw[0],raw[1])).all()); + for(int k=0; k<2; ++k) VERIFY(m2(k) == DenseIndex(raw[k])); + for(int k=0; k<2; ++k) VERIFY(a2(k) == DenseIndex(raw[k])); + } + { + Matrix m(raw), m2( (DenseIndex(raw[0])), (DenseIndex(raw[1])) ); + Array a(raw), a2( (DenseIndex(raw[0])), (DenseIndex(raw[1])) ); + for(int k=0; k<2; ++k) VERIFY(m(k) == raw[k]); + for(int k=0; k<2; ++k) VERIFY(a(k) == raw[k]); + VERIFY_IS_EQUAL(m,(Matrix(raw[0],raw[1]))); + VERIFY((a==Array(raw[0],raw[1])).all()); + for(int k=0; k<2; ++k) VERIFY(m2(k) == DenseIndex(raw[k])); + for(int k=0; k<2; ++k) VERIFY(a2(k) == DenseIndex(raw[k])); } { - Matrix m(raw); - Array a(raw); + Matrix m(raw), m1(raw[0]), m2( (DenseIndex(raw[0])) ); + Array a(raw), a1(raw[0]), a2( (DenseIndex(raw[0])) ); VERIFY(m(0) == raw[0]); VERIFY(a(0) == raw[0]); + VERIFY(m1(0) == raw[0]); + VERIFY(a1(0) == raw[0]); + VERIFY(m2(0) == DenseIndex(raw[0])); + VERIFY(a2(0) == DenseIndex(raw[0])); VERIFY_IS_EQUAL(m,(Matrix(raw[0]))); VERIFY((a==Array(raw[0])).all()); } @@ -233,9 +250,10 @@ void test_basicstuff() } CALL_SUBTEST_1(fixedSizeMatrixConstruction()); - CALL_SUBTEST_1(fixedSizeMatrixConstruction()); + CALL_SUBTEST_1(fixedSizeMatrixConstruction()); CALL_SUBTEST_1(fixedSizeMatrixConstruction()); CALL_SUBTEST_1(fixedSizeMatrixConstruction()); + CALL_SUBTEST_1(fixedSizeMatrixConstruction()); CALL_SUBTEST_1(fixedSizeMatrixConstruction()); CALL_SUBTEST_2(casting()); -- cgit v1.2.3 From 7f15f27a9e239bd386541b9ff33ace35fb2a6994 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Jul 2014 17:01:34 +0200 Subject: Workaround ambiguous call of init1 with MSVC. --- Eigen/src/Core/PlainObjectBase.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/PlainObjectBase.h b/Eigen/src/Core/PlainObjectBase.h index b968b325e..82d96eb92 100644 --- a/Eigen/src/Core/PlainObjectBase.h +++ b/Eigen/src/Core/PlainObjectBase.h @@ -721,7 +721,7 @@ class PlainObjectBase : public internal::dense_xpr_base::type EIGEN_STRONG_INLINE void _init1(const Index& val0, typename internal::enable_if< (!internal::is_same::value) && Base::SizeAtCompileTime==1 - && internal::is_convertible::value,T>::type* = 0) + && internal::is_convertible::value,T*>::type* = 0) { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 1) m_storage.data()[0] = Scalar(val0); -- cgit v1.2.3 From d1e9f39a9aac85eeec2e17f00aef24cb35537be3 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Jul 2014 18:28:19 +0200 Subject: Ambiguous call fixes for clang. --- Eigen/src/Core/PlainObjectBase.h | 11 ++++++++--- test/basicstuff.cpp | 10 ++++++++-- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Core/PlainObjectBase.h b/Eigen/src/Core/PlainObjectBase.h index 82d96eb92..f8292c793 100644 --- a/Eigen/src/Core/PlainObjectBase.h +++ b/Eigen/src/Core/PlainObjectBase.h @@ -693,7 +693,11 @@ class PlainObjectBase : public internal::dense_xpr_base::type template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE void _init2(const Index& val0, const Index& val1, typename internal::enable_if<(!internal::is_same::value) && Base::SizeAtCompileTime==2,T1>::type* = 0) + EIGEN_STRONG_INLINE void _init2(const Index& val0, const Index& val1, + typename internal::enable_if< (!internal::is_same::value) + && (internal::is_same::value) + && (internal::is_same::value) + && Base::SizeAtCompileTime==2,T1>::type* = 0) { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 2) m_storage.data()[0] = Scalar(val0); @@ -719,8 +723,9 @@ class PlainObjectBase : public internal::dense_xpr_base::type template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void _init1(const Index& val0, - typename internal::enable_if< (!internal::is_same::value) - && Base::SizeAtCompileTime==1 + typename internal::enable_if< (!internal::is_same::value) + && (internal::is_same::value) + && Base::SizeAtCompileTime==1 && internal::is_convertible::value,T*>::type* = 0) { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(PlainObjectBase, 1) diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 70cefe19f..3d7d74d4e 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -211,7 +211,10 @@ void fixedSizeMatrixConstruction() for(int k=0; k<2; ++k) VERIFY(a2(k) == DenseIndex(raw[k])); } { - Matrix m(raw), m2( (DenseIndex(raw[0])), (DenseIndex(raw[1])) ); + Matrix m(raw), + m2( (DenseIndex(raw[0])), (DenseIndex(raw[1])) ), + m3( (int(raw[0])), (int(raw[1])) ), + m4( (float(raw[0])), (float(raw[1])) ); Array a(raw), a2( (DenseIndex(raw[0])), (DenseIndex(raw[1])) ); for(int k=0; k<2; ++k) VERIFY(m(k) == raw[k]); for(int k=0; k<2; ++k) VERIFY(a(k) == raw[k]); @@ -219,9 +222,11 @@ void fixedSizeMatrixConstruction() VERIFY((a==Array(raw[0],raw[1])).all()); for(int k=0; k<2; ++k) VERIFY(m2(k) == DenseIndex(raw[k])); for(int k=0; k<2; ++k) VERIFY(a2(k) == DenseIndex(raw[k])); + for(int k=0; k<2; ++k) VERIFY(m3(k) == int(raw[k])); + for(int k=0; k<2; ++k) VERIFY(m4(k) == float(raw[k])); } { - Matrix m(raw), m1(raw[0]), m2( (DenseIndex(raw[0])) ); + Matrix m(raw), m1(raw[0]), m2( (DenseIndex(raw[0])) ), m3( (int(raw[0])) ); Array a(raw), a1(raw[0]), a2( (DenseIndex(raw[0])) ); VERIFY(m(0) == raw[0]); VERIFY(a(0) == raw[0]); @@ -229,6 +234,7 @@ void fixedSizeMatrixConstruction() VERIFY(a1(0) == raw[0]); VERIFY(m2(0) == DenseIndex(raw[0])); VERIFY(a2(0) == DenseIndex(raw[0])); + VERIFY(m3(0) == int(raw[0])); VERIFY_IS_EQUAL(m,(Matrix(raw[0]))); VERIFY((a==Array(raw[0])).all()); } -- cgit v1.2.3 From 2c625ec9ba64c6c09217190b494c9b08efb499e5 Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Tue, 22 Jul 2014 20:46:03 +0000 Subject: Simplification of some Altivec constants, reuse existing constants and avoid loading from RAM esp in the case of p16uc_COMPLEX_TRANSPOSE* --- Eigen/src/Core/arch/AltiVec/Complex.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h index 73fa62643..13b874d0c 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h @@ -16,13 +16,14 @@ namespace internal { static Packet4ui p4ui_CONJ_XOR = vec_mergeh((Packet4ui)p4i_ZERO, (Packet4ui)p4f_ZERO_);//{ 0x00000000, 0x80000000, 0x00000000, 0x80000000 }; static Packet16uc p16uc_COMPLEX_RE = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 0), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 2), 8);//{ 0,1,2,3, 0,1,2,3, 8,9,10,11, 8,9,10,11 }; -static Packet16uc p16uc_COMPLEX_IM = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 1), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 }; +static Packet16uc p16uc_COMPLEX_IM = vec_sld(p16uc_DUPLICATE, (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 }; static Packet16uc p16uc_COMPLEX_REV = vec_sld(p16uc_REVERSE, p16uc_REVERSE, 8);//{ 4,5,6,7, 0,1,2,3, 12,13,14,15, 8,9,10,11 }; static Packet16uc p16uc_COMPLEX_REV2 = vec_sld(p16uc_FORWARD, p16uc_FORWARD, 8);//{ 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 }; -static Packet16uc p16uc_PSET_HI = (Packet16uc) vec_mergeh((Packet4ui) vec_splat((Packet4ui)p16uc_FORWARD, 0), (Packet4ui) vec_splat((Packet4ui)p16uc_FORWARD, 1));//{ 0,1,2,3, 4,5,6,7, 0,1,2,3, 4,5,6,7 }; -static Packet16uc p16uc_PSET_LO = (Packet16uc) vec_mergeh((Packet4ui) vec_splat((Packet4ui)p16uc_FORWARD, 2), (Packet4ui) vec_splat((Packet4ui)p16uc_FORWARD, 3));//{ 8,9,10,11, 12,13,14,15, 8,9,10,11, 12,13,14,15 }; -static Packet16uc p16uc_COMPLEX_TRANSPOSE_0 = { 0,1,2,3, 4,5,6,7, 16,17,18,19, 20,21,22,23}; -static Packet16uc p16uc_COMPLEX_TRANSPOSE_1 = { 8,9,10,11, 12,13,14,15, 24,25,26,27, 28,29,30,31}; +static Packet16uc p16uc_PSET_HI = (Packet16uc) vec_mergeh((Packet4ui)p16uc_COMPLEX_RE, (Packet4ui)p16uc_COMPLEX_IM);//{ 0,1,2,3, 4,5,6,7, 0,1,2,3, 4,5,6,7 }; +static Packet16uc p16uc_PSET_LO = (Packet16uc) vec_mergel((Packet4ui)p16uc_COMPLEX_RE, (Packet4ui)p16uc_COMPLEX_IM);//{ 8,9,10,11, 12,13,14,15, 8,9,10,11, 12,13,14,15 }; +static Packet16uc p16uc_COMPLEX_MASK16 = vec_sld((Packet16uc)p4i_ZERO, vec_splat((Packet16uc) vec_abs(p4i_MINUS16), 3), 8);//{ 0,0,0,0, 0,0,0,0, 16,16,16,16, 16,16,16,16}; +static Packet16uc p16uc_COMPLEX_TRANSPOSE_0 = vec_add(p16uc_PSET_HI, p16uc_COMPLEX_MASK16);//{ 0,1,2,3, 4,5,6,7, 16,17,18,19, 20,21,22,23}; +static Packet16uc p16uc_COMPLEX_TRANSPOSE_1 = vec_add(p16uc_PSET_LO, p16uc_COMPLEX_MASK16);//{ 8,9,10,11, 12,13,14,15, 24,25,26,27, 28,29,30,31}; //---------- float ---------- struct Packet2cf -- cgit v1.2.3 From a0a87410d0cd76e930279f1a7c81c6da90148e59 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 24 Jul 2014 22:08:10 +0200 Subject: Fix bug #61: gemm was broken since we changed the blocking order --- blas/level3_impl.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/blas/level3_impl.h b/blas/level3_impl.h index 07dbc22ff..a05872666 100644 --- a/blas/level3_impl.h +++ b/blas/level3_impl.h @@ -56,7 +56,7 @@ int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScal else matrix(c, *m, *n, *ldc) *= beta; } - internal::gemm_blocking_space blocking(*m,*n,*k); + internal::gemm_blocking_space blocking(*m,*n,*k,true); int code = OP(*opa) | (OP(*opb) << 2); func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0); -- cgit v1.2.3 From 5f3d542b8a526d2fb8a179eb68f4a6d35f4f85a3 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Fri, 25 Jul 2014 13:34:03 +0100 Subject: Fix typo in MatrixExponential noticed by Markos. --- unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 05df5a102..160120d03 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -176,8 +176,8 @@ void matrix_exp_pade17(const MatrixType &A, MatrixType &U, MatrixType &V) const MatrixType A4 = A2 * A2; const MatrixType A6 = A4 * A2; const MatrixType A8 = A4 * A4; - V = b[17] * m_tmp1 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage - matrixType tmp = A8 * V; + V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage + MatrixType tmp = A8 * V; tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); U.noalias() = A * tmp; -- cgit v1.2.3 From ba694ce8cff79521850eee3285606589925880a5 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 30 Jul 2014 09:32:35 +0200 Subject: add missing delete operator overloads --- Eigen/src/Core/util/Memory.h | 2 ++ test/dynalloc.cpp | 31 +++++++++++++++++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 38990566d..810ee786b 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -703,6 +703,8 @@ template class aligned_stack_memory_handler } \ void operator delete(void * ptr) throw() { Eigen::internal::conditional_aligned_free(ptr); } \ void operator delete[](void * ptr) throw() { Eigen::internal::conditional_aligned_free(ptr); } \ + void operator delete(void * ptr, std::size_t /* sz */) throw() { Eigen::internal::conditional_aligned_free(ptr); } \ + void operator delete[](void * ptr, std::size_t /* sz */) throw() { Eigen::internal::conditional_aligned_free(ptr); } \ /* in-place new and delete. since (at least afaik) there is no actual */ \ /* memory allocated we can safely let the default implementation handle */ \ /* this particular case. */ \ diff --git a/test/dynalloc.cpp b/test/dynalloc.cpp index c98cc80f0..4f53ea6f0 100644 --- a/test/dynalloc.cpp +++ b/test/dynalloc.cpp @@ -93,6 +93,32 @@ template void check_dynaligned() } } +template void check_custom_new_delete() +{ + { + T* t = new T; + delete t; + } + + { + std::size_t N = internal::random(1,10); + T* t = new T[N]; + delete[] t; + } + +#ifdef EIGEN_ALIGN + { + T* t = static_cast((T::operator new)(sizeof(T))); + (T::operator delete)(t, sizeof(T)); + } + + { + T* t = static_cast((T::operator new)(sizeof(T))); + (T::operator delete)(t); + } +#endif +} + void test_dynalloc() { // low level dynamic memory allocation @@ -109,6 +135,11 @@ void test_dynalloc() CALL_SUBTEST(check_dynaligned() ); CALL_SUBTEST(check_dynaligned() ); CALL_SUBTEST(check_dynaligned() ); + + CALL_SUBTEST( check_custom_new_delete() ); + CALL_SUBTEST( check_custom_new_delete() ); + CALL_SUBTEST( check_custom_new_delete() ); + CALL_SUBTEST( check_custom_new_delete() ); } // check static allocation, who knows ? -- cgit v1.2.3 From d79516660c1f42ae0719ed353614d0977bd40153 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 31 Jul 2014 16:43:19 +0200 Subject: Make loadMarket use the sparse-matrix index type, thus enabling loading huge matrices. --- unsupported/Eigen/src/SparseExtra/MarketIO.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h index 1c40d3f7c..25ff4228d 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -133,6 +133,7 @@ template bool loadMarket(SparseMatrixType& mat, const std::string& filename) { typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::Index Index; std::ifstream input(filename.c_str(),std::ios::in); if(!input) return false; @@ -142,11 +143,11 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename) bool readsizes = false; - typedef Triplet T; + typedef Triplet T; std::vector elements; - int M(-1), N(-1), NNZ(-1); - int count = 0; + Index M(-1), N(-1), NNZ(-1); + Index count = 0; while(input.getline(buffer, maxBuffersize)) { // skip comments @@ -169,7 +170,7 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename) } else { - int i(-1), j(-1); + Index i(-1), j(-1); Scalar value; if( internal::GetMarketLine(line, M, N, i, j, value) ) { -- cgit v1.2.3 From db76193bc7dba3c58503bc45f0c12f23ecd296e2 Mon Sep 17 00:00:00 2001 From: Benjamin Chrétien Date: Fri, 1 Aug 2014 14:41:49 +0200 Subject: Fix typo in PermutationMatrix (doc). --- Eigen/src/Core/PermutationMatrix.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index 009873c78..8aa4c8bc5 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -262,7 +262,7 @@ class PermutationBase : public EigenBase * * \param SizeAtCompileTime the number of rows/cols, or Dynamic * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it. - * \param StorageIndexType the interger type of the indices + * \param StorageIndexType the integer type of the indices * * This class represents a permutation matrix, internally stored as a vector of integers. * -- cgit v1.2.3 From 6f58a41097986d53d2ece3525550a7f7657d05ed Mon Sep 17 00:00:00 2001 From: Benjamin Chrétien Date: Fri, 1 Aug 2014 15:35:45 +0200 Subject: Fix typos in Ref.h (doc). --- Eigen/src/Core/Ref.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Core/Ref.h b/Eigen/src/Core/Ref.h index cd6d949c4..12c9c09f4 100644 --- a/Eigen/src/Core/Ref.h +++ b/Eigen/src/Core/Ref.h @@ -39,10 +39,10 @@ template& x); * \endcode * - * In the in-out case, the input argument must satisfies the constraints of the actual Ref<> type, otherwise a compilation issue will be triggered. + * In the in-out case, the input argument must satisfy the constraints of the actual Ref<> type, otherwise a compilation issue will be triggered. * By default, a Ref can reference any dense vector expression of float having a contiguous memory layout. - * Likewise, a Ref can reference any column major dense matrix expression of float whose column's elements are contiguously stored with - * the possibility to have a constant space inbetween each column, i.e.: the inner stride mmust be equal to 1, but the outer-stride (or leading dimension), + * Likewise, a Ref can reference any column-major dense matrix expression of float whose column's elements are contiguously stored with + * the possibility to have a constant space in-between each column, i.e. the inner stride must be equal to 1, but the outer stride (or leading dimension) * can be greater than the number of rows. * * In the const case, if the input expression does not match the above requirement, then it is evaluated into a temporary before being passed to the function. -- cgit v1.2.3 From c53f88297c8dc0a1622258df701bdc864ff9ee76 Mon Sep 17 00:00:00 2001 From: Benjamin Chrétien Date: Fri, 1 Aug 2014 15:43:47 +0200 Subject: Fix more typos in Ref.h (doc). --- Eigen/src/Core/Ref.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Core/Ref.h b/Eigen/src/Core/Ref.h index 12c9c09f4..92614c6e2 100644 --- a/Eigen/src/Core/Ref.h +++ b/Eigen/src/Core/Ref.h @@ -19,17 +19,17 @@ template object can represent either a const expression or a l-value: * \code * // in-out argument: @@ -58,15 +58,15 @@ template > x); * foo3(A.row()); // OK * \endcode - * The downside here is that the function foo3 might be significantly slower than foo1 because it won't be able to exploit vectorization, and will involved more - * expensive address computations even if the input is contiguously stored in memory. To overcome this issue, one might propose to overloads internally calling a + * The downside here is that the function foo3 might be significantly slower than foo1 because it won't be able to exploit vectorization, and will involve more + * expensive address computations even if the input is contiguously stored in memory. To overcome this issue, one might propose to overload internally calling a * template function, e.g.: * \code * // in the .h: -- cgit v1.2.3 From 3e59163a24c5ff1a8009887cb5d5e091ae6df540 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 2 Aug 2014 02:47:30 +0200 Subject: Fix bug #850: workaround MSVC 2008 weird compilation bug --- Eigen/src/Core/PlainObjectBase.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/PlainObjectBase.h b/Eigen/src/Core/PlainObjectBase.h index f8292c793..69f34bd3e 100644 --- a/Eigen/src/Core/PlainObjectBase.h +++ b/Eigen/src/Core/PlainObjectBase.h @@ -708,7 +708,9 @@ class PlainObjectBase : public internal::dense_xpr_base::type EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void _init1(Index size, typename internal::enable_if::value,T>::type* = 0) { - EIGEN_STATIC_ASSERT(bool(NumTraits::IsInteger), + // NOTE MSVC 2008 complains if we directly put bool(NumTraits::IsInteger) as the EIGEN_STATIC_ASSERT argument. + const bool is_integer = NumTraits::IsInteger; + EIGEN_STATIC_ASSERT(is_integer, FLOATING_POINT_ARGUMENT_PASSED__INTEGER_WAS_EXPECTED) resize(size); } -- cgit v1.2.3 From 953ec08089ac7c5fc3c8f11236f3c546c88b5853 Mon Sep 17 00:00:00 2001 From: Kolja Brix Date: Sat, 2 Aug 2014 18:39:15 +0200 Subject: Correct GMRES: * Fix error in calculation of residual at restart. * Use relative residual as stopping criterion. * Improve documentation. --- unsupported/Eigen/src/IterativeSolvers/GMRES.h | 156 +++++++++++++------------ 1 file changed, 80 insertions(+), 76 deletions(-) diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index c8c84069e..67498705b 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -11,7 +11,7 @@ #ifndef EIGEN_GMRES_H #define EIGEN_GMRES_H -namespace Eigen { +namespace Eigen { namespace internal { @@ -27,11 +27,11 @@ namespace internal { * \param iters on input: maximum number of iterations to perform * on output: number of iterations performed * \param restart number of iterations for a restart - * \param tol_error on input: residual tolerance + * \param tol_error on input: relative residual tolerance * on output: residuum achieved * - * \sa IterativeMethods::bicgstab() - * + * \sa IterativeMethods::bicgstab() + * * * For references, please see: * @@ -70,18 +70,24 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition const int m = mat.rows(); - VectorType p0 = rhs - mat*x; + // residual and preconditioned residual + const VectorType p0 = rhs - mat*x; VectorType r0 = precond.solve(p0); - + + const RealScalar r0Norm = r0.norm(); + // is initial guess already good enough? - if(abs(r0.norm()) < tol) { - return true; + if(r0Norm == 0) { + tol_error=0; + return true; } + // storage for Hessenberg matrix and Householder data + FMatrixType H = FMatrixType::Zero(m, restart + 1); VectorType w = VectorType::Zero(restart + 1); - - FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix VectorType tau = VectorType::Zero(restart + 1); + + // storage for Jacobi rotations std::vector < JacobiRotation < Scalar > > G(restart); // generate first Householder vector @@ -112,11 +118,10 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition } if (v.tail(m - k).norm() != 0.0) { - if (k <= restart) { // generate new Householder vector - VectorType e(m - k - 1); + VectorType e(m - k - 1); RealScalar beta; v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta); H.col(k).tail(m - k - 1) = e; @@ -125,78 +130,77 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data()); } - } + } - if (k > 1) { - for (int i = 0; i < k - 1; ++i) { - // apply old Givens rotations to v - v.applyOnTheLeft(i, i + 1, G[i].adjoint()); - } - } + if (k > 1) { + for (int i = 0; i < k - 1; ++i) { + // apply old Givens rotations to v + v.applyOnTheLeft(i, i + 1, G[i].adjoint()); + } + } - if (k ().solveInPlace(y); + // solve upper triangular system + VectorType y = w.head(k); + H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y); - // use Horner-like scheme to calculate solution vector - VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); + // use Horner-like scheme to calculate solution vector + VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); - // apply Householder reflection H_{k} to x_new - x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); + // apply Householder reflection H_{k} to x_new + x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); - for (int i = k - 2; i >= 0; --i) { - x_new += y(i) * VectorType::Unit(m, i); - // apply Householder reflection H_{i} to x_new - x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); - } + for (int i = k - 2; i >= 0; --i) { + x_new += y(i) * VectorType::Unit(m, i); + // apply Householder reflection H_{i} to x_new + x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); + } - x += x_new; + x += x_new; - if (stop) { - return true; - } else { - k=0; + if (stop) { + return true; + } else { + k=0; - // reset data for a restart r0 = rhs - mat * x; - VectorType p0=mat*x; - VectorType p1=precond.solve(p0); - r0 = rhs - p1; -// r0_sqnorm = r0.squaredNorm(); - w = VectorType::Zero(restart + 1); - H = FMatrixType::Zero(m, restart + 1); - tau = VectorType::Zero(restart + 1); + // reset data for restart + const VectorType p0 = rhs - mat*x; + r0 = precond.solve(p0); - // generate first Householder vector - RealScalar beta; - r0.makeHouseholder(e, tau.coeffRef(0), beta); - w(0)=(Scalar) beta; - H.bottomLeftCorner(m - 1, 1) = e; + // clear Hessenberg matrix and Householder data + H = FMatrixType::Zero(m, restart + 1); + w = VectorType::Zero(restart + 1); + tau = VectorType::Zero(restart + 1); - } + // generate first Householder vector + RealScalar beta; + r0.makeHouseholder(e, tau.coeffRef(0), beta); + w(0)=(Scalar) beta; + H.bottomLeftCorner(m - 1, 1) = e; - } + } + } } - + return false; } @@ -230,7 +234,7 @@ struct traits > * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations * and NumTraits::epsilon() for the tolerance. - * + * * This class can be used as the direct solver classes. Here is a typical usage example: * \code * int n = 10000; @@ -244,7 +248,7 @@ struct traits > * // update b, and solve again * x = solver.solve(b); * \endcode - * + * * By default the iterations start with x=0 as an initial guess of the solution. * One can control the start using the solveWithGuess() method. Here is a step by * step execution example starting with a random guess and printing the evolution @@ -260,7 +264,7 @@ struct traits > * } while (solver.info()!=Success && i<100); * \endcode * Note that such a step by step excution is slightly slower. - * + * * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner */ template< typename _MatrixType, typename _Preconditioner> @@ -272,10 +276,10 @@ class GMRES : public IterativeSolverBase > using Base::m_iterations; using Base::m_info; using Base::m_isInitialized; - + private: int m_restart; - + public: typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; @@ -289,10 +293,10 @@ public: GMRES() : Base(), m_restart(30) {} /** Initialize the solver with matrix \a A for further \c Ax=b solving. - * + * * This constructor is a shortcut for the default constructor followed * by a call to compute(). - * + * * \warning this class stores a reference to the matrix A as well as some * precomputed values that depend on it. Therefore, if \a A is changed * this class becomes invalid. Call compute() to update it with the new @@ -301,16 +305,16 @@ public: GMRES(const MatrixType& A) : Base(A), m_restart(30) {} ~GMRES() {} - + /** Get the number of iterations after that a restart is performed. */ int get_restart() { return m_restart; } - + /** Set the number of iterations after that a restart is performed. * \param restart number of iterations for a restarti, default is 30. */ void set_restart(const int restart) { m_restart=restart; } - + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A * \a x0 as an initial solution. * @@ -326,17 +330,17 @@ public: return internal::solve_retval_with_guess (*this, b.derived(), x0); } - + /** \internal */ template void _solveWithGuess(const Rhs& b, Dest& x) const - { + { bool failed = false; for(int j=0; j Date: Mon, 4 Aug 2014 12:46:00 +0200 Subject: Memory allocated on the stack is freed at the function exit, so reduce iteration count to avoid stack overflow --- test/dynalloc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/dynalloc.cpp b/test/dynalloc.cpp index 4f53ea6f0..1190eb9cd 100644 --- a/test/dynalloc.cpp +++ b/test/dynalloc.cpp @@ -55,7 +55,7 @@ void check_aligned_new() void check_aligned_stack_alloc() { - for(int i = 1; i < 1000; i++) + for(int i = 1; i < 400; i++) { ei_declare_aligned_stack_constructed_variable(float,p,i,0); VERIFY(size_t(p)%ALIGNMENT==0); -- cgit v1.2.3 From 50085d2c28bcd4e3b8e1eb2e11dec60312cb8ac3 Mon Sep 17 00:00:00 2001 From: Silvio Traversaro Date: Thu, 7 Aug 2014 15:48:53 +0000 Subject: FindEigen3.cmake: Add reading hints of Eigen directory location from environment variables EIGEN3_ROOT and EIGEN3_ROOT_DIR . --- cmake/FindEigen3.cmake | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/cmake/FindEigen3.cmake b/cmake/FindEigen3.cmake index 9c546a05d..cea1afeab 100644 --- a/cmake/FindEigen3.cmake +++ b/cmake/FindEigen3.cmake @@ -9,6 +9,12 @@ # EIGEN3_FOUND - system has eigen lib with correct version # EIGEN3_INCLUDE_DIR - the eigen include directory # EIGEN3_VERSION - eigen version +# +# This module reads hints about search locations from +# the following enviroment variables: +# +# EIGEN3_ROOT +# EIGEN3_ROOT_DIR # Copyright (c) 2006, 2007 Montel Laurent, # Copyright (c) 2008, 2009 Gael Guennebaud, @@ -62,6 +68,9 @@ if (EIGEN3_INCLUDE_DIR) else (EIGEN3_INCLUDE_DIR) find_path(EIGEN3_INCLUDE_DIR NAMES signature_of_eigen3_matrix_library + HINTS + ENV EIGEN3_ROOT + ENV EIGEN3_ROOT_DIR PATHS ${CMAKE_INSTALL_PREFIX}/include ${KDE4_INCLUDE_DIR} -- cgit v1.2.3 From 6a3423da800234645ea3a192507a0ee2fe317d6c Mon Sep 17 00:00:00 2001 From: Vladimir Chalupecky Date: Wed, 20 Aug 2014 05:42:22 +0000 Subject: Fix uninitialized variable warning in SparseQR --- Eigen/src/SparseQR/SparseQR.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 4c6553bf2..66dccfc43 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -447,7 +447,7 @@ void SparseQR::factorize(const MatrixType& mat) } } // End update current column - Scalar tau; + Scalar tau = RealScalar(0); RealScalar beta = 0; if(nonzeroCol < diagSize) @@ -461,7 +461,6 @@ void SparseQR::factorize(const MatrixType& mat) for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) { - tau = RealScalar(0); beta = numext::real(c0); tval(Qidx(0)) = 1; } -- cgit v1.2.3 From eeadc06e838a737492625cc3e4d5d5555bb40ff7 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Wed, 20 Aug 2014 16:39:25 +0200 Subject: EIGEN_EXCEPTIONS was not defined in test/main.h, therefore all VERIFY_RAISES_ASSERT tests were not enabled --- Eigen/Core | 2 +- test/main.h | 6 +++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/Eigen/Core b/Eigen/Core index 9a73fe37b..776b7faf3 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -42,7 +42,7 @@ #define EIGEN_USING_STD_MATH(FUNC) using std::FUNC; #endif -#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(__CUDA_ARCH__) +#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(__CUDA_ARCH__) && !defined(EIGEN_EXCEPTIONS) #define EIGEN_EXCEPTIONS #endif diff --git a/test/main.h b/test/main.h index 3295dcb71..376232cf2 100644 --- a/test/main.h +++ b/test/main.h @@ -76,6 +76,10 @@ namespace Eigen #define EIGEN_DEFAULT_IO_FORMAT IOFormat(4, 0, " ", "\n", "", "", "", "") +#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(__CUDA_ARCH__) + #define EIGEN_EXCEPTIONS +#endif + #ifndef EIGEN_NO_ASSERTION_CHECKING namespace Eigen @@ -172,7 +176,7 @@ namespace Eigen #ifndef VERIFY_RAISES_ASSERT #define VERIFY_RAISES_ASSERT(a) \ - std::cout << "Can't VERIFY_RAISES_ASSERT( " #a " ) with exceptions disabled"; + std::cout << "Can't VERIFY_RAISES_ASSERT( " #a " ) with exceptions disabled\n"; #endif #if !defined(__CUDACC__) -- cgit v1.2.3 From 9c0aa81fbfc4abd0dc297d75305b9d579dad2754 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 21 Aug 2014 10:49:09 +0200 Subject: bug #854: fix numerical issue in SelfAdjointEigenSolver::computeDirect for 3x3 matrices. The tolerance to detect stable cross products was too optimistic. Add respective unit tests. --- Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 18 ++++++++--- test/eigensolver_selfadjoint.cpp | 43 ++++++++++++++++++-------- 2 files changed, 44 insertions(+), 17 deletions(-) diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index fc8ecaa6f..a6bbdac6b 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -605,7 +605,6 @@ template struct direct_selfadjoint_eigenvalues::epsilon(); - safeNorm2 *= safeNorm2; if((eivals(2)-eivals(0))<=Eigen::NumTraits::epsilon()) { eivecs.setIdentity(); @@ -619,7 +618,7 @@ template struct direct_selfadjoint_eigenvalues d1 ? 2 : 0; - d0 = d0 > d1 ? d1 : d0; + d0 = d0 > d1 ? d0 : d1; tmp.diagonal().array () -= eivals(k); VectorType cross; @@ -627,19 +626,25 @@ template struct direct_selfadjoint_eigenvaluessafeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { // the input matrix and/or the eigenvaues probably contains some inf/NaN, @@ -659,12 +664,16 @@ template struct direct_selfadjoint_eigenvalues::epsilon()) + { eivecs.col(1) = eivecs.col(k).unitOrthogonal(); + } else { - n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm(); + n = (cross = eivecs.col(k).cross(tmp.row(0))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(1) = cross / sqrt(n); + } else { n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm(); @@ -678,13 +687,14 @@ template struct direct_selfadjoint_eigenvalues void selfadjointeigensolver(const MatrixType& m) MatrixType a = MatrixType::Random(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols); MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; + MatrixType symmC = symmA; + + // randomly nullify some rows/columns + { + Index count = 1;//internal::random(-cols,cols); + for(Index k=0; k(0,cols-1); + symmA.row(i).setZero(); + symmA.col(i).setZero(); + } + } + symmA.template triangularView().setZero(); + symmC.template triangularView().setZero(); MatrixType b = MatrixType::Random(rows,cols); MatrixType b1 = MatrixType::Random(rows,cols); @@ -40,7 +54,7 @@ template void selfadjointeigensolver(const MatrixType& m) SelfAdjointEigenSolver eiDirect; eiDirect.computeDirect(symmA); // generalized eigen pb - GeneralizedSelfAdjointEigenSolver eiSymmGen(symmA, symmB); + GeneralizedSelfAdjointEigenSolver eiSymmGen(symmC, symmB); VERIFY_IS_EQUAL(eiSymm.info(), Success); VERIFY((symmA.template selfadjointView() * eiSymm.eigenvectors()).isApprox( @@ -57,27 +71,28 @@ template void selfadjointeigensolver(const MatrixType& m) VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues()); // generalized eigen problem Ax = lBx - eiSymmGen.compute(symmA, symmB,Ax_lBx); + eiSymmGen.compute(symmC, symmB,Ax_lBx); VERIFY_IS_EQUAL(eiSymmGen.info(), Success); - VERIFY((symmA.template selfadjointView() * eiSymmGen.eigenvectors()).isApprox( + VERIFY((symmC.template selfadjointView() * eiSymmGen.eigenvectors()).isApprox( symmB.template selfadjointView() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); // generalized eigen problem BAx = lx - eiSymmGen.compute(symmA, symmB,BAx_lx); + eiSymmGen.compute(symmC, symmB,BAx_lx); VERIFY_IS_EQUAL(eiSymmGen.info(), Success); - VERIFY((symmB.template selfadjointView() * (symmA.template selfadjointView() * eiSymmGen.eigenvectors())).isApprox( + VERIFY((symmB.template selfadjointView() * (symmC.template selfadjointView() * eiSymmGen.eigenvectors())).isApprox( (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); // generalized eigen problem ABx = lx - eiSymmGen.compute(symmA, symmB,ABx_lx); + eiSymmGen.compute(symmC, symmB,ABx_lx); VERIFY_IS_EQUAL(eiSymmGen.info(), Success); - VERIFY((symmA.template selfadjointView() * (symmB.template selfadjointView() * eiSymmGen.eigenvectors())).isApprox( + VERIFY((symmC.template selfadjointView() * (symmB.template selfadjointView() * eiSymmGen.eigenvectors())).isApprox( (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); + eiSymm.compute(symmC); MatrixType sqrtSymmA = eiSymm.operatorSqrt(); - VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView()), sqrtSymmA*sqrtSymmA); - VERIFY_IS_APPROX(sqrtSymmA, symmA.template selfadjointView()*eiSymm.operatorInverseSqrt()); + VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView()), sqrtSymmA*sqrtSymmA); + VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView()*eiSymm.operatorInverseSqrt()); MatrixType id = MatrixType::Identity(rows, cols); VERIFY_IS_APPROX(id.template selfadjointView().operatorNorm(), RealScalar(1)); @@ -95,9 +110,9 @@ template void selfadjointeigensolver(const MatrixType& m) VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt()); // test Tridiagonalization's methods - Tridiagonalization tridiag(symmA); + Tridiagonalization tridiag(symmC); // FIXME tridiag.matrixQ().adjoint() does not work - VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); + VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); // Test computation of eigenvalues from tridiagonal matrix if(rows > 1) @@ -111,8 +126,8 @@ template void selfadjointeigensolver(const MatrixType& m) if (rows > 1) { // Test matrix with NaN - symmA(0,0) = std::numeric_limits::quiet_NaN(); - SelfAdjointEigenSolver eiSymmNaN(symmA); + symmC(0,0) = std::numeric_limits::quiet_NaN(); + SelfAdjointEigenSolver eiSymmNaN(symmC); VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence); } } @@ -122,8 +137,10 @@ void test_eigensolver_selfadjoint() int s = 0; for(int i = 0; i < g_repeat; i++) { // very important to test 3x3 and 2x2 matrices since we provide special paths for them + CALL_SUBTEST_1( selfadjointeigensolver(Matrix2f()) ); CALL_SUBTEST_1( selfadjointeigensolver(Matrix2d()) ); CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) ); + CALL_SUBTEST_1( selfadjointeigensolver(Matrix3d()) ); CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) ); s = internal::random(1,EIGEN_TEST_MAX_SIZE/4); CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) ); -- cgit v1.2.3 From b49ef99617b11404007db283ce68f5dd5b2eff0e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 21 Aug 2014 22:14:25 +0200 Subject: Do not apply the preconditioner before starting the iterations as this might destroy a very good initial guess. --- Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index dc524c225..27824b9d5 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -39,7 +39,6 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, int maxIters = iters; int n = mat.cols(); - x = precond.solve(x); VectorType r = rhs - mat * x; VectorType r0 = r; -- cgit v1.2.3 From 2e50289ba323740ef89eb51a74e00e645a8eb9be Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 26 Aug 2014 12:54:19 +0200 Subject: bug #861: enable posix_memalign with PGI --- Eigen/src/Core/util/Memory.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 810ee786b..30133ba67 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -64,7 +64,7 @@ // Currently, let's include it only on unix systems: #if defined(__unix__) || defined(__unix) #include - #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0) + #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || (defined __PGI) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0) #define EIGEN_HAS_POSIX_MEMALIGN 1 #endif #endif -- cgit v1.2.3 From be3477e2069cbfe7158a00c11e265d81d7e88e7e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 26 Aug 2014 12:52:29 +0200 Subject: bug #857: workaround MSVC compilation issue. --- Eigen/src/SparseCore/SparseDenseProduct.h | 9 --------- Eigen/src/SparseCore/SparseMatrixBase.h | 3 ++- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index 4a7813296..d40e966c1 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h @@ -318,15 +318,6 @@ class DenseTimeSparseProduct DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&); }; -// sparse * dense -template -template -inline const typename SparseDenseProductReturnType::Type -SparseMatrixBase::operator*(const MatrixBase &other) const -{ - return typename SparseDenseProductReturnType::Type(derived(), other.derived()); -} - } // end namespace Eigen #endif // EIGEN_SPARSEDENSEPRODUCT_H diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index 1050cf3f1..fb5025049 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -358,7 +358,8 @@ template class SparseMatrixBase : public EigenBase /** sparse * dense (returns a dense object unless it is an outer product) */ template const typename SparseDenseProductReturnType::Type - operator*(const MatrixBase &other) const; + operator*(const MatrixBase &other) const + { return typename SparseDenseProductReturnType::Type(derived(), other.derived()); } /** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */ SparseSymmetricPermutationProduct twistedBy(const PermutationMatrix& perm) const -- cgit v1.2.3 From 25a3e65a685f66f04b83003f12cd97b91e646c20 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 26 Aug 2014 23:32:32 +0200 Subject: In SparseQR, calling factorize() without analyzePattern() was broken. --- Eigen/src/SparseQR/SparseQR.h | 16 ++++++++++++++-- test/sparseqr.cpp | 2 ++ 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 66dccfc43..002b4824b 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -75,7 +75,7 @@ class SparseQR typedef Matrix ScalarVector; typedef PermutationMatrix PermutationType; public: - SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) + SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { } /** Construct a QR factorization of the matrix \a mat. @@ -84,7 +84,7 @@ class SparseQR * * \sa compute() */ - SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) + SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { compute(mat); } @@ -262,6 +262,7 @@ class SparseQR IndexVector m_etree; // Column elimination tree IndexVector m_firstRowElt; // First element in each row bool m_isQSorted; // whether Q is sorted or not + bool m_isEtreeOk; // whether the elimination tree match the initial input matrix template friend struct SparseQR_QProduct; template friend struct SparseQRMatrixQReturnType; @@ -297,6 +298,7 @@ void SparseQR::analyzePattern(const MatrixType& mat) // Compute the column elimination tree of the permuted matrix m_outputPerm_c = m_perm_c.inverse(); internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + m_isEtreeOk = true; m_R.resize(m, n); m_Q.resize(m, diagSize); @@ -330,6 +332,15 @@ void SparseQR::factorize(const MatrixType& mat) Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q ScalarVector tval(m); // The dense vector used to compute the current column RealScalar pivotThreshold = m_threshold; + + m_R.setZero(); + m_Q.setZero(); + if(!m_isEtreeOk) + { + m_outputPerm_c = m_perm_c.inverse(); + internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + m_isEtreeOk = true; + } m_pmat = mat; m_pmat.uncompress(); // To have the innerNonZeroPtr allocated @@ -513,6 +524,7 @@ void SparseQR::factorize(const MatrixType& mat) // Recompute the column elimination tree internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data()); + m_isEtreeOk = false; } } diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp index 1fe4a98ee..8e6887dd3 100644 --- a/test/sparseqr.cpp +++ b/test/sparseqr.cpp @@ -54,6 +54,8 @@ template void test_sparseqr_scalar() b = dA * DenseVector::Random(A.cols()); solver.compute(A); + if(internal::random(0,1)>0.5) + solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change. if (solver.info() != Success) { std::cerr << "sparse QR factorization failed\n"; -- cgit v1.2.3 From 0ba490cf80d9c389de410beaa3551b2a2a72a801 Mon Sep 17 00:00:00 2001 From: Georg Drenkhahn Date: Fri, 22 Aug 2014 12:13:07 +0200 Subject: Fixed CMakeLists.txt files to prevent CMake 3.0.0 warnings about deprecated LOCATION target property. Small whitespace cleanup in CMakelLists.txt. --- doc/examples/CMakeLists.txt | 4 +--- doc/snippets/CMakeLists.txt | 6 ++---- doc/special_examples/CMakeLists.txt | 16 +++++----------- unsupported/doc/examples/CMakeLists.txt | 4 +--- unsupported/doc/snippets/CMakeLists.txt | 4 +--- 5 files changed, 10 insertions(+), 24 deletions(-) diff --git a/doc/examples/CMakeLists.txt b/doc/examples/CMakeLists.txt index 71b272a15..08cf8efd7 100644 --- a/doc/examples/CMakeLists.txt +++ b/doc/examples/CMakeLists.txt @@ -6,12 +6,10 @@ foreach(example_src ${examples_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - get_target_property(example_executable - ${example} LOCATION) add_custom_command( TARGET ${example} POST_BUILD - COMMAND ${example_executable} + COMMAND ${example} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out ) add_dependencies(all_examples ${example}) diff --git a/doc/snippets/CMakeLists.txt b/doc/snippets/CMakeLists.txt index 92a22ea61..1135900cf 100644 --- a/doc/snippets/CMakeLists.txt +++ b/doc/snippets/CMakeLists.txt @@ -14,12 +14,10 @@ foreach(snippet_src ${snippets_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - get_target_property(compile_snippet_executable - ${compile_snippet_target} LOCATION) add_custom_command( TARGET ${compile_snippet_target} POST_BUILD - COMMAND ${compile_snippet_executable} + COMMAND ${compile_snippet_target} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out ) add_dependencies(all_snippets ${compile_snippet_target}) @@ -27,4 +25,4 @@ foreach(snippet_src ${snippets_SRCS}) PROPERTIES OBJECT_DEPENDS ${snippet_src}) endforeach(snippet_src) -ei_add_target_property(compile_tut_arithmetic_transpose_aliasing COMPILE_FLAGS -DEIGEN_NO_DEBUG) \ No newline at end of file +ei_add_target_property(compile_tut_arithmetic_transpose_aliasing COMPILE_FLAGS -DEIGEN_NO_DEBUG) diff --git a/doc/special_examples/CMakeLists.txt b/doc/special_examples/CMakeLists.txt index 45e2339e3..aab80a55d 100644 --- a/doc/special_examples/CMakeLists.txt +++ b/doc/special_examples/CMakeLists.txt @@ -1,4 +1,3 @@ - if(NOT EIGEN_TEST_NOQT) find_package(Qt4) if(QT4_FOUND) @@ -6,17 +5,16 @@ if(NOT EIGEN_TEST_NOQT) endif() endif(NOT EIGEN_TEST_NOQT) - if(QT4_FOUND) add_executable(Tutorial_sparse_example Tutorial_sparse_example.cpp Tutorial_sparse_example_details.cpp) target_link_libraries(Tutorial_sparse_example ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO} ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY}) - + add_custom_command( TARGET Tutorial_sparse_example POST_BUILD COMMAND Tutorial_sparse_example ARGS ${CMAKE_CURRENT_BINARY_DIR}/../html/Tutorial_sparse_example.jpeg ) - + add_dependencies(all_examples Tutorial_sparse_example) endif(QT4_FOUND) @@ -26,15 +24,11 @@ if(EIGEN_COMPILER_SUPPORT_CPP11) target_link_libraries(random_cpp11 ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) add_dependencies(all_examples random_cpp11) ei_add_target_property(random_cpp11 COMPILE_FLAGS "-std=c++11") - - get_target_property(random_cpp11_exec - random_cpp11 LOCATION) + add_custom_command( TARGET random_cpp11 POST_BUILD - COMMAND ${random_cpp11_exec} + COMMAND random_cpp11 ARGS >${CMAKE_CURRENT_BINARY_DIR}/random_cpp11.out ) - - -endif() \ No newline at end of file +endif() diff --git a/unsupported/doc/examples/CMakeLists.txt b/unsupported/doc/examples/CMakeLists.txt index 978f9afd8..c47646dfc 100644 --- a/unsupported/doc/examples/CMakeLists.txt +++ b/unsupported/doc/examples/CMakeLists.txt @@ -10,12 +10,10 @@ FOREACH(example_src ${examples_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(example_${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - GET_TARGET_PROPERTY(example_executable - example_${example} LOCATION) ADD_CUSTOM_COMMAND( TARGET example_${example} POST_BUILD - COMMAND ${example_executable} + COMMAND example_${example} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out ) ADD_DEPENDENCIES(unsupported_examples example_${example}) diff --git a/unsupported/doc/snippets/CMakeLists.txt b/unsupported/doc/snippets/CMakeLists.txt index 4a4157933..f0c5cc2a8 100644 --- a/unsupported/doc/snippets/CMakeLists.txt +++ b/unsupported/doc/snippets/CMakeLists.txt @@ -14,12 +14,10 @@ FOREACH(snippet_src ${snippets_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - GET_TARGET_PROPERTY(compile_snippet_executable - ${compile_snippet_target} LOCATION) ADD_CUSTOM_COMMAND( TARGET ${compile_snippet_target} POST_BUILD - COMMAND ${compile_snippet_executable} + COMMAND ${compile_snippet_target} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out ) ADD_DEPENDENCIES(unsupported_snippets ${compile_snippet_target}) -- cgit v1.2.3 From c3e408047427a12669720b64397e080956786829 Mon Sep 17 00:00:00 2001 From: Freddie Witherden Date: Wed, 27 Aug 2014 15:24:51 +0100 Subject: Allow LevenbergMarquardt to work with non-standard types. --- .../Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h | 8 +++++--- .../Eigen/src/NonLinearOptimization/LevenbergMarquardt.h | 12 +++++++----- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h index 51dd1d3c4..7cebe4e06 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h @@ -144,11 +144,13 @@ class LevenbergMarquardt : internal::no_assignment_operator /** Sets the default parameters */ void resetParameters() - { + { + using std::sqrt; + m_factor = 100.; m_maxfev = 400; - m_ftol = std::sqrt(NumTraits::epsilon()); - m_xtol = std::sqrt(NumTraits::epsilon()); + m_ftol = sqrt(NumTraits::epsilon()); + m_xtol = sqrt(NumTraits::epsilon()); m_gtol = 0. ; m_epsfcn = 0. ; } diff --git a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h index bfeb26fc9..ecb8dccf4 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h @@ -55,8 +55,8 @@ public: Parameters() : factor(Scalar(100.)) , maxfev(400) - , ftol(std::sqrt(NumTraits::epsilon())) - , xtol(std::sqrt(NumTraits::epsilon())) + , ftol(sqrt_(NumTraits::epsilon())) + , xtol(sqrt_(NumTraits::epsilon())) , gtol(Scalar(0.)) , epsfcn(Scalar(0.)) {} Scalar factor; @@ -72,7 +72,7 @@ public: LevenbergMarquardtSpace::Status lmder1( FVectorType &x, - const Scalar tol = std::sqrt(NumTraits::epsilon()) + const Scalar tol = sqrt_(NumTraits::epsilon()) ); LevenbergMarquardtSpace::Status minimize(FVectorType &x); @@ -83,12 +83,12 @@ public: FunctorType &functor, FVectorType &x, Index *nfev, - const Scalar tol = std::sqrt(NumTraits::epsilon()) + const Scalar tol = sqrt_(NumTraits::epsilon()) ); LevenbergMarquardtSpace::Status lmstr1( FVectorType &x, - const Scalar tol = std::sqrt(NumTraits::epsilon()) + const Scalar tol = sqrt_(NumTraits::epsilon()) ); LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x); @@ -109,6 +109,8 @@ public: Scalar lm_param(void) { return par; } private: + static Scalar sqrt_(const Scalar& x) { using std::sqrt; return sqrt(x); } + FunctorType &functor; Index n; Index m; -- cgit v1.2.3 From 1ed9e2d0044a5c35052965ea53e4f905279a06f1 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 29 Aug 2014 11:55:03 +0200 Subject: In sparse matrix product, enable sorted insertion when doing two transposition is defenitely not optimal. --- .../SparseCore/ConservativeSparseSparseProduct.h | 82 ++++++++++++---------- 1 file changed, 45 insertions(+), 37 deletions(-) diff --git a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h index 5c320e2d2..608044a95 100644 --- a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h +++ b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h @@ -15,7 +15,7 @@ namespace Eigen { namespace internal { template -static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) +static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res, bool sortedInsertion = false) { typedef typename remove_all::type::Scalar Scalar; typedef typename remove_all::type::Index Index; @@ -64,53 +64,51 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r values[i] += x * y; } } - - // unordered insertion - for(Index k=0; k use a quick sort - // otherwise => loop through the entire vector - // In order to avoid to perform an expensive log2 when the - // result is clearly very sparse we use a linear bound up to 200. - //if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); + // unordered insertion for(Index k=0; k use a quick sort + // otherwise => loop through the entire vector + // In order to avoid to perform an expensive log2 when the + // result is clearly very sparse we use a linear bound up to 200. + if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); + for(Index k=0; k RowMajorMatrix; typedef SparseMatrix ColMajorMatrix; ColMajorMatrix resCol(lhs.rows(),rhs.cols()); - internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); - // sort the non zeros: - RowMajorMatrix resRow(resCol); - res = resRow; + // FIXME, the following heuristic is probably not very good. + if(lhs.rows()>=rhs.cols()) + { + // perform sorted insertion + internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol, true); + res = resCol; + } + else + { + // ressort to transpose to sort the entries + internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol, false); + RowMajorMatrix resRow(resCol); + res = resRow; + } } }; -- cgit v1.2.3 From 460662cbcc89c378b9ea097220b77d0eea9551ff Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 29 Aug 2014 14:18:23 +0200 Subject: Fix SparseVector::coeffRef(i,j) and add missing SparseVector::insert*Unordered --- Eigen/src/SparseCore/SparseVector.h | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/Eigen/src/SparseCore/SparseVector.h b/Eigen/src/SparseCore/SparseVector.h index 7e15c814b..0b1b389ce 100644 --- a/Eigen/src/SparseCore/SparseVector.h +++ b/Eigen/src/SparseCore/SparseVector.h @@ -109,7 +109,7 @@ class SparseVector inline Scalar& coeffRef(Index row, Index col) { eigen_assert(IsColVector ? (col==0 && row>=0 && row=0 && col Date: Fri, 29 Aug 2014 14:19:03 +0200 Subject: Optimization in sparse-sparse matrix products for small ones --- Eigen/src/SparseCore/ConservativeSparseSparseProduct.h | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h index 608044a95..8067565f9 100644 --- a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h +++ b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h @@ -24,10 +24,10 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r Index rows = lhs.innerSize(); Index cols = rhs.outerSize(); eigen_assert(lhs.outerSize() == rhs.innerSize()); - - std::vector mask(rows,false); - Matrix values(rows); - Matrix indices(rows); + + ei_declare_aligned_stack_constructed_variable(bool, mask, rows, 0); + ei_declare_aligned_stack_constructed_variable(Scalar, values, rows, 0); + ei_declare_aligned_stack_constructed_variable(Index, indices, rows, 0); // estimate the number of non zero entries // given a rhs column containing Y non zeros, we assume that the respective Y columns @@ -77,7 +77,7 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r else { // alternative ordered insertion code: - const Index t200 = rows/(log2(200)*1.39); + const Index t200 = rows/11; // 11 == (log2(200)*1.39) const Index t = (rows*100)/139; // FIXME reserve nnz non zeros @@ -88,7 +88,7 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r // result is clearly very sparse we use a linear bound up to 200. if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); + if(nnz>1) std::sort(indices,indices+nnz); for(Index k=0; k RowMajorMatrix; - typedef SparseMatrix ColMajorMatrix; + typedef SparseMatrix ColMajorMatrixAux; + typedef typename sparse_eval::type ColMajorMatrix; + ColMajorMatrix resCol(lhs.rows(),rhs.cols()); // FIXME, the following heuristic is probably not very good. if(lhs.rows()>=rhs.cols()) { // perform sorted insertion internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol, true); - res = resCol; + res.swap(resCol); } else { -- cgit v1.2.3