From 47913da94a27e98a9115f85b2a530b6c14a10b8f Mon Sep 17 00:00:00 2001 From: Jose Luis Date: Fri, 2 May 2014 00:22:51 +0200 Subject: [PATCH] Upgrade embedded Eigen to 3.2.1 --- otherlibs/eigen3/Eigen/Cholesky | 3 +- otherlibs/eigen3/Eigen/CholmodSupport | 3 +- otherlibs/eigen3/Eigen/Core | 4 +- otherlibs/eigen3/Eigen/Eigen2Support | 18 +- otherlibs/eigen3/Eigen/Eigenvalues | 3 +- otherlibs/eigen3/Eigen/Geometry | 3 +- otherlibs/eigen3/Eigen/Householder | 3 +- otherlibs/eigen3/Eigen/IterativeLinearSolvers | 3 +- otherlibs/eigen3/Eigen/Jacobi | 3 +- otherlibs/eigen3/Eigen/LU | 3 +- otherlibs/eigen3/Eigen/LeastSquares | 3 +- otherlibs/eigen3/Eigen/MetisSupport | 3 +- otherlibs/eigen3/Eigen/OrderingMethods | 3 +- otherlibs/eigen3/Eigen/PaStiXSupport | 3 +- otherlibs/eigen3/Eigen/PardisoSupport | 3 +- otherlibs/eigen3/Eigen/QR | 3 +- otherlibs/eigen3/Eigen/SPQRSupport | 3 +- otherlibs/eigen3/Eigen/SVD | 3 +- otherlibs/eigen3/Eigen/Sparse | 3 +- otherlibs/eigen3/Eigen/SparseCholesky | 3 +- otherlibs/eigen3/Eigen/SparseCore | 3 +- otherlibs/eigen3/Eigen/SparseLU | 3 +- otherlibs/eigen3/Eigen/SparseQR | 3 +- otherlibs/eigen3/Eigen/SuperLUSupport | 3 +- otherlibs/eigen3/Eigen/UmfPackSupport | 3 +- otherlibs/eigen3/Eigen/src/Cholesky/LDLT.h | 52 ++-- .../Eigen/src/CholmodSupport/CholmodSupport.h | 23 +- otherlibs/eigen3/Eigen/src/Core/Array.h | 5 +- .../eigen3/Eigen/src/Core/BooleanRedux.h | 8 +- otherlibs/eigen3/Eigen/src/Core/EigenBase.h | 30 --- .../eigen3/Eigen/src/Core/GeneralProduct.h | 17 +- otherlibs/eigen3/Eigen/src/Core/IO.h | 9 +- otherlibs/eigen3/Eigen/src/Core/Matrix.h | 5 +- otherlibs/eigen3/Eigen/src/Core/MatrixBase.h | 45 ++++ .../eigen3/Eigen/src/Core/PermutationMatrix.h | 3 +- .../eigen3/Eigen/src/Core/PlainObjectBase.h | 16 +- otherlibs/eigen3/Eigen/src/Core/Ref.h | 5 +- otherlibs/eigen3/Eigen/src/Core/StableNorm.h | 27 +- otherlibs/eigen3/Eigen/src/Core/Transpose.h | 4 +- .../eigen3/Eigen/src/Core/VectorwiseOp.h | 5 +- .../Eigen/src/Core/arch/SSE/MathFunctions.h | 11 + .../Eigen/src/Core/arch/SSE/PacketMath.h | 11 +- .../Core/products/GeneralBlockPanelKernel.h | 6 + .../src/Core/products/GeneralMatrixVector.h | 15 +- .../Core/products/SelfadjointMatrixVector.h | 6 +- .../eigen3/Eigen/src/Core/util/Constants.h | 8 +- otherlibs/eigen3/Eigen/src/Core/util/Macros.h | 9 +- otherlibs/eigen3/Eigen/src/Core/util/Memory.h | 13 +- .../eigen3/Eigen/src/Eigen2Support/SVD.h | 3 +- .../eigen3/Eigen/src/Geometry/AlignedBox.h | 3 +- .../eigen3/Eigen/src/Geometry/EulerAngles.h | 2 +- .../eigen3/Eigen/src/Geometry/Quaternion.h | 27 +- .../eigen3/Eigen/src/Geometry/Transform.h | 4 +- .../eigen3/Eigen/src/QR/ColPivHouseholderQR.h | 2 +- .../Eigen/src/QR/FullPivHouseholderQR.h | 53 ++-- .../src/SPQRSupport/SuiteSparseQRSupport.h | 18 +- otherlibs/eigen3/Eigen/src/SVD/JacobiSVD.h | 24 +- .../ConservativeSparseSparseProduct.h | 34 +-- .../Eigen/src/SparseCore/MappedSparseMatrix.h | 2 + .../eigen3/Eigen/src/SparseCore/SparseBlock.h | 4 + .../Eigen/src/SparseCore/SparseColEtree.h | 6 +- .../Eigen/src/SparseCore/SparseDenseProduct.h | 4 +- .../Eigen/src/SparseCore/SparseMatrix.h | 21 +- .../Eigen/src/SparseCore/SparseMatrixBase.h | 4 +- .../Eigen/src/SparseCore/SparsePermutation.h | 4 +- .../Eigen/src/SparseCore/SparseProduct.h | 5 +- .../SparseSparseProductWithPruning.h | 13 +- .../eigen3/Eigen/src/SparseCore/SparseUtil.h | 2 +- .../eigen3/Eigen/src/SparseLU/SparseLU.h | 8 +- .../Eigen/src/SparseLU/SparseLU_Memory.h | 43 ++-- .../eigen3/Eigen/src/SparseQR/SparseQR.h | 13 +- .../eigen3/Eigen/src/plugins/BlockMethods.h | 240 ++++++++++-------- .../Eigen/src/plugins/MatrixCwiseBinaryOps.h | 4 +- .../eigen3/unsupported/Eigen/AdolcForward | 3 +- .../eigen3/unsupported/Eigen/AlignedVector3 | 3 +- .../eigen3/unsupported/Eigen/ArpackSupport | 3 +- otherlibs/eigen3/unsupported/Eigen/AutoDiff | 3 +- otherlibs/eigen3/unsupported/Eigen/BVH | 3 +- otherlibs/eigen3/unsupported/Eigen/FFT | 3 +- .../eigen3/unsupported/Eigen/IterativeSolvers | 3 +- .../eigen3/unsupported/Eigen/KroneckerProduct | 3 +- .../unsupported/Eigen/LevenbergMarquardt | 3 +- .../eigen3/unsupported/Eigen/MPRealSupport | 3 +- .../eigen3/unsupported/Eigen/MatrixFunctions | 3 +- .../unsupported/Eigen/MoreVectorization | 3 +- .../unsupported/Eigen/NonLinearOptimization | 3 +- .../eigen3/unsupported/Eigen/NumericalDiff | 3 +- .../eigen3/unsupported/Eigen/OpenGLSupport | 10 +- .../eigen3/unsupported/Eigen/Polynomials | 3 +- otherlibs/eigen3/unsupported/Eigen/SVD | 3 +- otherlibs/eigen3/unsupported/Eigen/Skyline | 3 +- .../eigen3/unsupported/Eigen/SparseExtra | 3 +- otherlibs/eigen3/unsupported/Eigen/Splines | 3 +- .../IterativeSolvers/ConstrainedConjGrad.h | 10 +- otherlibs/eigen3/unsupported/doc/Overview.dox | 3 +- otherlibs/eigen3/unsupported/test/FFTW.cpp | 13 +- 96 files changed, 567 insertions(+), 493 deletions(-) diff --git a/otherlibs/eigen3/Eigen/Cholesky b/otherlibs/eigen3/Eigen/Cholesky index ccbf6fe72d..f727f5d89c 100644 --- a/otherlibs/eigen3/Eigen/Cholesky +++ b/otherlibs/eigen3/Eigen/Cholesky @@ -5,8 +5,7 @@ #include "src/Core/util/DisableStupidWarnings.h" -/** \defgroup Cholesky_Module Cholesky module - * \ingroup eigen_grp +/** \defgroup Cholesky_Module Cholesky module * * * diff --git a/otherlibs/eigen3/Eigen/CholmodSupport b/otherlibs/eigen3/Eigen/CholmodSupport index b83c594e52..745b884e74 100644 --- a/otherlibs/eigen3/Eigen/CholmodSupport +++ b/otherlibs/eigen3/Eigen/CholmodSupport @@ -10,8 +10,7 @@ extern "C" { } /** \ingroup Support_modules - * \defgroup CholmodSupport_Module CholmodSupport module - * \ingroup eigen_grp + * \defgroup CholmodSupport_Module CholmodSupport module * * This module provides an interface to the Cholmod library which is part of the suitesparse package. * It provides the two following main factorization classes: diff --git a/otherlibs/eigen3/Eigen/Core b/otherlibs/eigen3/Eigen/Core index fb0df47e4d..9131cc3fc9 100644 --- a/otherlibs/eigen3/Eigen/Core +++ b/otherlibs/eigen3/Eigen/Core @@ -238,9 +238,7 @@ using std::size_t; // gcc 4.6.0 wants std:: for ptrdiff_t using std::ptrdiff_t; -/** \defgroup eigen_grp The Eigen3 library - \defgroup Core_Module Core module - * \ingroup eigen_grp +/** \defgroup Core_Module Core module * This is the main module of Eigen providing dense matrix and vector support * (both fixed and dynamic size) with all the features corresponding to a BLAS library * and much more... diff --git a/otherlibs/eigen3/Eigen/Eigen2Support b/otherlibs/eigen3/Eigen/Eigen2Support index fb33ffa1a5..6aa009d209 100644 --- a/otherlibs/eigen3/Eigen/Eigen2Support +++ b/otherlibs/eigen3/Eigen/Eigen2Support @@ -14,13 +14,25 @@ #error Eigen2 support must be enabled by defining EIGEN2_SUPPORT before including any Eigen header #endif +#ifndef EIGEN_NO_EIGEN2_DEPRECATED_WARNING + +#if defined(__GNUC__) || defined(__INTEL_COMPILER) || defined(__clang__) +#warning "Eigen2 support is deprecated in Eigen 3.2.x and it will be removed in Eigen 3.3. (Define EIGEN_NO_EIGEN2_DEPRECATED_WARNING to disable this warning)" +#else +#pragma message ("Eigen2 support is deprecated in Eigen 3.2.x and it will be removed in Eigen 3.3. (Define EIGEN_NO_EIGEN2_DEPRECATED_WARNING to disable this warning)") +#endif + +#endif // EIGEN_NO_EIGEN2_DEPRECATED_WARNING + #include "src/Core/util/DisableStupidWarnings.h" /** \ingroup Support_modules - * \defgroup Eigen2Support_Module Eigen2 support module - * \ingroup eigen_grp - * This module provides a couple of deprecated functions improving the compatibility with Eigen2. + * \defgroup Eigen2Support_Module Eigen2 support module + * + * \warning Eigen2 support is deprecated in Eigen 3.2.x and it will be removed in Eigen 3.3. * + * This module provides a couple of deprecated functions improving the compatibility with Eigen2. + * * To use it, define EIGEN2_SUPPORT before including any Eigen header * \code * #define EIGEN2_SUPPORT diff --git a/otherlibs/eigen3/Eigen/Eigenvalues b/otherlibs/eigen3/Eigen/Eigenvalues index fbdc9d3fb1..53c5a73a27 100644 --- a/otherlibs/eigen3/Eigen/Eigenvalues +++ b/otherlibs/eigen3/Eigen/Eigenvalues @@ -11,8 +11,7 @@ #include "LU" #include "Geometry" -/** \defgroup Eigenvalues_Module Eigenvalues module - * \ingroup eigen_grp +/** \defgroup Eigenvalues_Module Eigenvalues module * * * diff --git a/otherlibs/eigen3/Eigen/Geometry b/otherlibs/eigen3/Eigen/Geometry index d458501ffc..efd9d4504c 100644 --- a/otherlibs/eigen3/Eigen/Geometry +++ b/otherlibs/eigen3/Eigen/Geometry @@ -13,8 +13,7 @@ #define M_PI 3.14159265358979323846 #endif -/** \defgroup Geometry_Module Geometry module - * \ingroup eigen_grp +/** \defgroup Geometry_Module Geometry module * * * diff --git a/otherlibs/eigen3/Eigen/Householder b/otherlibs/eigen3/Eigen/Householder index 9266bc43fb..6e348db5c4 100644 --- a/otherlibs/eigen3/Eigen/Householder +++ b/otherlibs/eigen3/Eigen/Householder @@ -5,8 +5,7 @@ #include "src/Core/util/DisableStupidWarnings.h" -/** \defgroup Householder_Module Householder module - * \ingroup eigen_grp +/** \defgroup Householder_Module Householder module * This module provides Householder transformations. * * \code diff --git a/otherlibs/eigen3/Eigen/IterativeLinearSolvers b/otherlibs/eigen3/Eigen/IterativeLinearSolvers index 3b50496b3c..0f4159dc19 100644 --- a/otherlibs/eigen3/Eigen/IterativeLinearSolvers +++ b/otherlibs/eigen3/Eigen/IterativeLinearSolvers @@ -7,8 +7,7 @@ #include "src/Core/util/DisableStupidWarnings.h" /** - * \defgroup IterativeLinearSolvers_Module IterativeLinearSolvers module - * \ingroup eigen_grp + * \defgroup IterativeLinearSolvers_Module IterativeLinearSolvers module * * This module currently provides iterative methods to solve problems of the form \c A \c x = \c b, where \c A is a squared matrix, usually very large and sparse. * Those solvers are accessible via the following classes: diff --git a/otherlibs/eigen3/Eigen/Jacobi b/otherlibs/eigen3/Eigen/Jacobi index 4031dfea12..ba8a4dc36a 100644 --- a/otherlibs/eigen3/Eigen/Jacobi +++ b/otherlibs/eigen3/Eigen/Jacobi @@ -5,8 +5,7 @@ #include "src/Core/util/DisableStupidWarnings.h" -/** \defgroup Jacobi_Module Jacobi module - * \ingroup eigen_grp +/** \defgroup Jacobi_Module Jacobi module * This module provides Jacobi and Givens rotations. * * \code diff --git a/otherlibs/eigen3/Eigen/LU b/otherlibs/eigen3/Eigen/LU index a8ed48abb0..db57955044 100644 --- a/otherlibs/eigen3/Eigen/LU +++ b/otherlibs/eigen3/Eigen/LU @@ -5,8 +5,7 @@ #include "src/Core/util/DisableStupidWarnings.h" -/** \defgroup LU_Module LU module - * \ingroup eigen_grp +/** \defgroup LU_Module LU module * This module includes %LU decomposition and related notions such as matrix inversion and determinant. * This module defines the following MatrixBase methods: * - MatrixBase::inverse() diff --git a/otherlibs/eigen3/Eigen/LeastSquares b/otherlibs/eigen3/Eigen/LeastSquares index 300f2eba9e..35137c25db 100644 --- a/otherlibs/eigen3/Eigen/LeastSquares +++ b/otherlibs/eigen3/Eigen/LeastSquares @@ -15,8 +15,7 @@ #include "Eigenvalues" #include "Geometry" -/** \defgroup LeastSquares_Module LeastSquares module - * \ingroup eigen_grp +/** \defgroup LeastSquares_Module LeastSquares module * This module provides linear regression and related features. * * \code diff --git a/otherlibs/eigen3/Eigen/MetisSupport b/otherlibs/eigen3/Eigen/MetisSupport index e29829f11f..6a113f7a87 100644 --- a/otherlibs/eigen3/Eigen/MetisSupport +++ b/otherlibs/eigen3/Eigen/MetisSupport @@ -11,8 +11,7 @@ extern "C" { /** \ingroup Support_modules - * \defgroup MetisSupport_Module MetisSupport module - * \ingroup eigen_grp + * \defgroup MetisSupport_Module MetisSupport module * * \code * #include diff --git a/otherlibs/eigen3/Eigen/OrderingMethods b/otherlibs/eigen3/Eigen/OrderingMethods index 0af901e625..7c0f1fffff 100644 --- a/otherlibs/eigen3/Eigen/OrderingMethods +++ b/otherlibs/eigen3/Eigen/OrderingMethods @@ -6,8 +6,7 @@ #include "src/Core/util/DisableStupidWarnings.h" /** - * \defgroup OrderingMethods_Module OrderingMethods module - * \ingroup eigen_grp + * \defgroup OrderingMethods_Module OrderingMethods module * * This module is currently for internal use only * diff --git a/otherlibs/eigen3/Eigen/PaStiXSupport b/otherlibs/eigen3/Eigen/PaStiXSupport index 455bf22e7e..7c616ee5ea 100644 --- a/otherlibs/eigen3/Eigen/PaStiXSupport +++ b/otherlibs/eigen3/Eigen/PaStiXSupport @@ -16,8 +16,7 @@ extern "C" { #endif /** \ingroup Support_modules - * \defgroup PaStiXSupport_Module PaStiXSupport module - * \ingroup eigen_grp + * \defgroup PaStiXSupport_Module PaStiXSupport module * * This module provides an interface to the PaSTiX library. * PaSTiX is a general \b supernodal, \b parallel and \b opensource sparse solver. diff --git a/otherlibs/eigen3/Eigen/PardisoSupport b/otherlibs/eigen3/Eigen/PardisoSupport index a5a586e353..99330ce7a7 100644 --- a/otherlibs/eigen3/Eigen/PardisoSupport +++ b/otherlibs/eigen3/Eigen/PardisoSupport @@ -10,8 +10,7 @@ #include /** \ingroup Support_modules - * \defgroup PardisoSupport_Module PardisoSupport module - * \ingroup eigen_grp + * \defgroup PardisoSupport_Module PardisoSupport module * * This module brings support for the Intel(R) MKL PARDISO direct sparse solvers. * diff --git a/otherlibs/eigen3/Eigen/QR b/otherlibs/eigen3/Eigen/QR index df21c36a09..ac5b026935 100644 --- a/otherlibs/eigen3/Eigen/QR +++ b/otherlibs/eigen3/Eigen/QR @@ -9,8 +9,7 @@ #include "Jacobi" #include "Householder" -/** \defgroup QR_Module QR module - * \ingroup eigen_grp +/** \defgroup QR_Module QR module * * * diff --git a/otherlibs/eigen3/Eigen/SPQRSupport b/otherlibs/eigen3/Eigen/SPQRSupport index c1cc0c6529..77016442ee 100644 --- a/otherlibs/eigen3/Eigen/SPQRSupport +++ b/otherlibs/eigen3/Eigen/SPQRSupport @@ -8,8 +8,7 @@ #include "SuiteSparseQR.hpp" /** \ingroup Support_modules - * \defgroup SPQRSupport_Module SuiteSparseQR module - * \ingroup eigen_grp + * \defgroup SPQRSupport_Module SuiteSparseQR module * * This module provides an interface to the SPQR library, which is part of the suitesparse package. * diff --git a/otherlibs/eigen3/Eigen/SVD b/otherlibs/eigen3/Eigen/SVD index 96d832c999..fd310017ad 100644 --- a/otherlibs/eigen3/Eigen/SVD +++ b/otherlibs/eigen3/Eigen/SVD @@ -7,8 +7,7 @@ #include "src/Core/util/DisableStupidWarnings.h" -/** \defgroup SVD_Module SVD module - * \ingroup eigen_grp +/** \defgroup SVD_Module SVD module * * * diff --git a/otherlibs/eigen3/Eigen/Sparse b/otherlibs/eigen3/Eigen/Sparse index 7865605b70..7cc9c09133 100644 --- a/otherlibs/eigen3/Eigen/Sparse +++ b/otherlibs/eigen3/Eigen/Sparse @@ -1,8 +1,7 @@ #ifndef EIGEN_SPARSE_MODULE_H #define EIGEN_SPARSE_MODULE_H -/** \defgroup Sparse_Module Sparse meta-module - * \ingroup eigen_grp +/** \defgroup Sparse_Module Sparse meta-module * * Meta-module including all related modules: * - \ref SparseCore_Module diff --git a/otherlibs/eigen3/Eigen/SparseCholesky b/otherlibs/eigen3/Eigen/SparseCholesky index 7220d17bf8..9f5056aa1a 100644 --- a/otherlibs/eigen3/Eigen/SparseCholesky +++ b/otherlibs/eigen3/Eigen/SparseCholesky @@ -16,8 +16,7 @@ #include "src/Core/util/DisableStupidWarnings.h" /** - * \defgroup SparseCholesky_Module SparseCholesky module - * \ingroup eigen_grp + * \defgroup SparseCholesky_Module SparseCholesky module * * This module currently provides two variants of the direct sparse Cholesky decomposition for selfadjoint (hermitian) matrices. * Those decompositions are accessible via the following classes: diff --git a/otherlibs/eigen3/Eigen/SparseCore b/otherlibs/eigen3/Eigen/SparseCore index a2d4b2fea7..9b5be5e15a 100644 --- a/otherlibs/eigen3/Eigen/SparseCore +++ b/otherlibs/eigen3/Eigen/SparseCore @@ -12,8 +12,7 @@ #include /** - * \defgroup SparseCore_Module SparseCore module - * \ingroup eigen_grp + * \defgroup SparseCore_Module SparseCore module * * This module provides a sparse matrix representation, and basic associatd matrix manipulations * and operations. diff --git a/otherlibs/eigen3/Eigen/SparseLU b/otherlibs/eigen3/Eigen/SparseLU index df117ca492..8527a49bd8 100644 --- a/otherlibs/eigen3/Eigen/SparseLU +++ b/otherlibs/eigen3/Eigen/SparseLU @@ -14,8 +14,7 @@ #include "SparseCore" /** - * \defgroup SparseLU_Module SparseLU module - * \ingroup eigen_grp + * \defgroup SparseLU_Module SparseLU module * This module defines a supernodal factorization of general sparse matrices. * The code is fully optimized for supernode-panel updates with specialized kernels. * Please, see the documentation of the SparseLU class for more details. diff --git a/otherlibs/eigen3/Eigen/SparseQR b/otherlibs/eigen3/Eigen/SparseQR index 62e46b9dca..4ee42065ee 100644 --- a/otherlibs/eigen3/Eigen/SparseQR +++ b/otherlibs/eigen3/Eigen/SparseQR @@ -5,8 +5,7 @@ #include "OrderingMethods" #include "src/Core/util/DisableStupidWarnings.h" -/** \defgroup SparseQR_Module SparseQR module - * \ingroup eigen_grp +/** \defgroup SparseQR_Module SparseQR module * \brief Provides QR decomposition for sparse matrices * * This module provides a simplicial version of the left-looking Sparse QR decomposition. diff --git a/otherlibs/eigen3/Eigen/SuperLUSupport b/otherlibs/eigen3/Eigen/SuperLUSupport index 81a225b8a4..575e14fbc2 100644 --- a/otherlibs/eigen3/Eigen/SuperLUSupport +++ b/otherlibs/eigen3/Eigen/SuperLUSupport @@ -29,8 +29,7 @@ typedef int int_t; namespace Eigen { struct SluMatrix; } /** \ingroup Support_modules - * \defgroup SuperLUSupport_Module SuperLUSupport module - * \ingroup eigen_grp + * \defgroup SuperLUSupport_Module SuperLUSupport module * * This module provides an interface to the SuperLU library. * It provides the following factorization class: diff --git a/otherlibs/eigen3/Eigen/UmfPackSupport b/otherlibs/eigen3/Eigen/UmfPackSupport index 4511b8bae6..984f64a841 100644 --- a/otherlibs/eigen3/Eigen/UmfPackSupport +++ b/otherlibs/eigen3/Eigen/UmfPackSupport @@ -10,8 +10,7 @@ extern "C" { } /** \ingroup Support_modules - * \defgroup UmfPackSupport_Module UmfPackSupport module - * \ingroup eigen_grp + * \defgroup UmfPackSupport_Module UmfPackSupport module * * This module provides an interface to the UmfPack library which is part of the suitesparse package. * It provides the following factorization class: diff --git a/otherlibs/eigen3/Eigen/src/Cholesky/LDLT.h b/otherlibs/eigen3/Eigen/src/Cholesky/LDLT.h index d19cb3968d..d026418f8a 100644 --- a/otherlibs/eigen3/Eigen/src/Cholesky/LDLT.h +++ b/otherlibs/eigen3/Eigen/src/Cholesky/LDLT.h @@ -16,7 +16,10 @@ namespace Eigen { namespace internal { -template struct LDLT_Traits; + template struct LDLT_Traits; + + // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef + enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite }; } /** \ingroup Cholesky_Module @@ -69,7 +72,12 @@ template class LDLT * The default constructor is useful in cases in which the user intends to * perform decompositions via LDLT::compute(const MatrixType&). */ - LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {} + LDLT() + : m_matrix(), + m_transpositions(), + m_sign(internal::ZeroSign), + m_isInitialized(false) + {} /** \brief Default Constructor with memory preallocation * @@ -81,6 +89,7 @@ template class LDLT : m_matrix(size, size), m_transpositions(size), m_temporary(size), + m_sign(internal::ZeroSign), m_isInitialized(false) {} @@ -93,6 +102,7 @@ template class LDLT : m_matrix(matrix.rows(), matrix.cols()), m_transpositions(matrix.rows()), m_temporary(matrix.rows()), + m_sign(internal::ZeroSign), m_isInitialized(false) { compute(matrix); @@ -139,7 +149,7 @@ template class LDLT inline bool isPositive() const { eigen_assert(m_isInitialized && "LDLT is not initialized."); - return m_sign == 1; + return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign; } #ifdef EIGEN2_SUPPORT @@ -153,7 +163,7 @@ template class LDLT inline bool isNegative(void) const { eigen_assert(m_isInitialized && "LDLT is not initialized."); - return m_sign == -1; + return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign; } /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -235,7 +245,7 @@ template class LDLT MatrixType m_matrix; TranspositionType m_transpositions; TmpMatrixType m_temporary; - int m_sign; + internal::SignMatrix m_sign; bool m_isInitialized; }; @@ -246,7 +256,7 @@ template struct ldlt_inplace; template<> struct ldlt_inplace { template - static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0) + static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) { using std::abs; typedef typename MatrixType::Scalar Scalar; @@ -258,8 +268,9 @@ template<> struct ldlt_inplace if (size <= 1) { transpositions.setIdentity(); - if(sign) - *sign = numext::real(mat.coeff(0,0))>0 ? 1:-1; + if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef; + else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef; + else sign = ZeroSign; return true; } @@ -284,7 +295,6 @@ template<> struct ldlt_inplace if(biggest_in_corner < cutoff) { for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; - if(sign) *sign = 0; break; } @@ -325,15 +335,15 @@ template<> struct ldlt_inplace } if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) A21 /= mat.coeffRef(k,k); - - if(sign) - { - // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right - int newSign = numext::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; - if(k == 0) - *sign = newSign; - else if(*sign != newSign) - *sign = 0; + + RealScalar realAkk = numext::real(mat.coeffRef(k,k)); + if (sign == PositiveSemiDef) { + if (realAkk < 0) sign = Indefinite; + } else if (sign == NegativeSemiDef) { + if (realAkk > 0) sign = Indefinite; + } else if (sign == ZeroSign) { + if (realAkk > 0) sign = PositiveSemiDef; + else if (realAkk < 0) sign = NegativeSemiDef; } } @@ -399,7 +409,7 @@ template<> struct ldlt_inplace template<> struct ldlt_inplace { template - static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0) + static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) { Transpose matt(mat); return ldlt_inplace::unblocked(matt, transpositions, temp, sign); @@ -445,7 +455,7 @@ LDLT& LDLT::compute(const MatrixType& a) m_isInitialized = false; m_temporary.resize(size); - internal::ldlt_inplace::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign); + internal::ldlt_inplace::unblocked(m_matrix, m_transpositions, m_temporary, m_sign); m_isInitialized = true; return *this; @@ -473,7 +483,7 @@ LDLT& LDLT::rankUpdate(const MatrixBase=0 ? 1 : -1; + m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef; m_isInitialized = true; } diff --git a/otherlibs/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h b/otherlibs/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h index 783324b0b2..c449960de4 100644 --- a/otherlibs/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/otherlibs/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h @@ -58,10 +58,12 @@ cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat) res.p = mat.outerIndexPtr(); res.i = mat.innerIndexPtr(); res.x = mat.valuePtr(); + res.z = 0; res.sorted = 1; if(mat.isCompressed()) { res.packed = 1; + res.nz = 0; } else { @@ -170,6 +172,7 @@ class CholmodBase : internal::noncopyable CholmodBase() : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) { + m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0); cholmod_start(&m_cholmod); } @@ -241,7 +244,7 @@ class CholmodBase : internal::noncopyable return internal::sparse_solve_retval(*this, b.derived()); } - /** Performs a symbolic decomposition on the sparcity of \a matrix. + /** Performs a symbolic decomposition on the sparsity pattern of \a matrix. * * This function is particularly useful when solving for several problems having the same structure. * @@ -265,7 +268,7 @@ class CholmodBase : internal::noncopyable /** Performs a numeric decomposition of \a matrix * - * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. + * The given matrix must have the same sparsity pattern as the matrix on which the symbolic decomposition has been performed. * * \sa analyzePattern() */ @@ -302,7 +305,7 @@ class CholmodBase : internal::noncopyable { this->m_info = NumericalIssue; } - // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.) + // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) dest = Matrix::Map(reinterpret_cast(x_cd->x),b.rows(),b.cols()); cholmod_free_dense(&x_cd, &m_cholmod); } @@ -323,7 +326,7 @@ class CholmodBase : internal::noncopyable { this->m_info = NumericalIssue; } - // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.) + // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) dest = viewAsEigen(*x_cs); cholmod_free_sparse(&x_cs, &m_cholmod); } @@ -365,8 +368,8 @@ class CholmodBase : internal::noncopyable * * This class allows to solve for A.X = B sparse linear problems via a simplicial LL^T Cholesky factorization * using the Cholmod library. - * This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Thefore, it has little practical interest. - * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Therefore, it has little practical interest. + * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices * X and B can be either dense or sparse. * * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> @@ -412,8 +415,8 @@ class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimpl * * This class allows to solve for A.X = B sparse linear problems via a simplicial LDL^T Cholesky factorization * using the Cholmod library. - * This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Thefore, it has little practical interest. - * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Therefore, it has little practical interest. + * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices * X and B can be either dense or sparse. * * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> @@ -458,7 +461,7 @@ class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimp * This class allows to solve for A.X = B sparse linear problems via a supernodal LL^T Cholesky factorization * using the Cholmod library. * This supernodal variant performs best on dense enough problems, e.g., 3D FEM, or very high order 2D FEM. - * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices * X and B can be either dense or sparse. * * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> @@ -501,7 +504,7 @@ class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSuper * \brief A general Cholesky factorization and solver based on Cholmod * * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization - * using the Cholmod library. The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * using the Cholmod library. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices * X and B can be either dense or sparse. * * This variant permits to change the underlying Cholesky method at runtime. diff --git a/otherlibs/eigen3/Eigen/src/Core/Array.h b/otherlibs/eigen3/Eigen/src/Core/Array.h index d6ffe743a1..0ab03eff0f 100644 --- a/otherlibs/eigen3/Eigen/src/Core/Array.h +++ b/otherlibs/eigen3/Eigen/src/Core/Array.h @@ -210,7 +210,7 @@ class Array : Base(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols()) { Base::_check_template_params(); - Base::resize(other.rows(), other.cols()); + Base::_resize_to_match(other); *this = other; } @@ -234,8 +234,7 @@ class Array friend struct internal::matrix_swap_impl; }; -/** \defgroup arraytypedefs Global array typedefs - * \ingroup eigen_grp +/** \defgroup arraytypedefs Global array typedefs * \ingroup Core_Module * * Eigen defines several typedef shortcuts for most common 1D and 2D array types. diff --git a/otherlibs/eigen3/Eigen/src/Core/BooleanRedux.h b/otherlibs/eigen3/Eigen/src/Core/BooleanRedux.h index 6e37e031a8..be9f48a8c7 100644 --- a/otherlibs/eigen3/Eigen/src/Core/BooleanRedux.h +++ b/otherlibs/eigen3/Eigen/src/Core/BooleanRedux.h @@ -29,9 +29,9 @@ struct all_unroller }; template -struct all_unroller +struct all_unroller { - static inline bool run(const Derived &mat) { return mat.coeff(0, 0); } + static inline bool run(const Derived &/*mat*/) { return true; } }; template @@ -55,9 +55,9 @@ struct any_unroller }; template -struct any_unroller +struct any_unroller { - static inline bool run(const Derived &mat) { return mat.coeff(0, 0); } + static inline bool run(const Derived & /*mat*/) { return false; } }; template diff --git a/otherlibs/eigen3/Eigen/src/Core/EigenBase.h b/otherlibs/eigen3/Eigen/src/Core/EigenBase.h index 2b8dd1b706..fadb45852f 100644 --- a/otherlibs/eigen3/Eigen/src/Core/EigenBase.h +++ b/otherlibs/eigen3/Eigen/src/Core/EigenBase.h @@ -126,36 +126,6 @@ Derived& DenseBase::operator-=(const EigenBase &other) return derived(); } -/** replaces \c *this by \c *this * \a other. - * - * \returns a reference to \c *this - */ -template -template -inline Derived& -MatrixBase::operator*=(const EigenBase &other) -{ - other.derived().applyThisOnTheRight(derived()); - return derived(); -} - -/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=(). - */ -template -template -inline void MatrixBase::applyOnTheRight(const EigenBase &other) -{ - other.derived().applyThisOnTheRight(derived()); -} - -/** replaces \c *this by \c *this * \a other. */ -template -template -inline void MatrixBase::applyOnTheLeft(const EigenBase &other) -{ - other.derived().applyThisOnTheLeft(derived()); -} - } // end namespace Eigen #endif // EIGEN_EIGENBASE_H diff --git a/otherlibs/eigen3/Eigen/src/Core/GeneralProduct.h b/otherlibs/eigen3/Eigen/src/Core/GeneralProduct.h index 578f96dffa..2a59d94645 100644 --- a/otherlibs/eigen3/Eigen/src/Core/GeneralProduct.h +++ b/otherlibs/eigen3/Eigen/src/Core/GeneralProduct.h @@ -11,7 +11,7 @@ #ifndef EIGEN_GENERAL_PRODUCT_H #define EIGEN_GENERAL_PRODUCT_H -namespace Eigen { +namespace Eigen { /** \class GeneralProduct * \ingroup Core_Module @@ -258,7 +258,7 @@ class GeneralProduct : public ProductBase, Lhs, Rhs> { template struct IsRowMajor : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {}; - + public: EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) @@ -267,7 +267,7 @@ class GeneralProduct EIGEN_STATIC_ASSERT((internal::is_same::value), YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) } - + struct set { template void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() = src; } }; struct add { template void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } }; struct sub { template void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } }; @@ -278,12 +278,12 @@ class GeneralProduct dst.const_cast_derived() += m_scale * src; } }; - + template inline void evalTo(Dest& dest) const { internal::outer_product_selector_run(*this, dest, set(), IsRowMajor()); } - + template inline void addTo(Dest& dest) const { internal::outer_product_selector_run(*this, dest, add(), IsRowMajor()); @@ -437,12 +437,12 @@ template<> struct gemv_selector bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0)); bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible; - + RhsScalar compatibleAlpha = get_factor::run(actualAlpha); ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(), evalToDest ? dest.data() : static_dest.data()); - + if(!evalToDest) { #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN @@ -587,7 +587,8 @@ MatrixBase::operator*(const MatrixBase &other) const EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), - INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) + INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) + EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) #ifdef EIGEN_DEBUG_PRODUCT internal::product_type::debug(); #endif diff --git a/otherlibs/eigen3/Eigen/src/Core/IO.h b/otherlibs/eigen3/Eigen/src/Core/IO.h index c8d5f6379b..8d4bc59e9d 100644 --- a/otherlibs/eigen3/Eigen/src/Core/IO.h +++ b/otherlibs/eigen3/Eigen/src/Core/IO.h @@ -185,21 +185,22 @@ std::ostream & print_matrix(std::ostream & s, const Derived& _m, const IOFormat& explicit_precision = fmt.precision; } + std::streamsize old_precision = 0; + if(explicit_precision) old_precision = s.precision(explicit_precision); + bool align_cols = !(fmt.flags & DontAlignCols); if(align_cols) { // compute the largest width - for(Index j = 1; j < m.cols(); ++j) + for(Index j = 0; j < m.cols(); ++j) for(Index i = 0; i < m.rows(); ++i) { std::stringstream sstr; - if(explicit_precision) sstr.precision(explicit_precision); + sstr.copyfmt(s); sstr << m.coeff(i,j); width = std::max(width, Index(sstr.str().length())); } } - std::streamsize old_precision = 0; - if(explicit_precision) old_precision = s.precision(explicit_precision); s << fmt.matPrefix; for(Index i = 0; i < m.rows(); ++i) { diff --git a/otherlibs/eigen3/Eigen/src/Core/Matrix.h b/otherlibs/eigen3/Eigen/src/Core/Matrix.h index bad8aa08b2..d7d0b5b9a4 100644 --- a/otherlibs/eigen3/Eigen/src/Core/Matrix.h +++ b/otherlibs/eigen3/Eigen/src/Core/Matrix.h @@ -304,7 +304,7 @@ class Matrix : Base(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols()) { Base::_check_template_params(); - Base::resize(other.rows(), other.cols()); + Base::_resize_to_match(other); // FIXME/CHECK: isn't *this = other.derived() more efficient. it allows to // go for pure _set() implementations, right? *this = other; @@ -347,8 +347,7 @@ class Matrix using Base::m_storage; }; -/** \defgroup matrixtypedefs Global matrix typedefs - * \ingroup eigen_grp +/** \defgroup matrixtypedefs Global matrix typedefs * * \ingroup Core_Module * diff --git a/otherlibs/eigen3/Eigen/src/Core/MatrixBase.h b/otherlibs/eigen3/Eigen/src/Core/MatrixBase.h index 9193b6abb3..344b38f2fc 100644 --- a/otherlibs/eigen3/Eigen/src/Core/MatrixBase.h +++ b/otherlibs/eigen3/Eigen/src/Core/MatrixBase.h @@ -510,6 +510,51 @@ template class MatrixBase {EIGEN_STATIC_ASSERT(std::ptrdiff_t(sizeof(typename OtherDerived::Scalar))==-1,YOU_CANNOT_MIX_ARRAYS_AND_MATRICES); return *this;} }; + +/*************************************************************************** +* Implementation of matrix base methods +***************************************************************************/ + +/** replaces \c *this by \c *this * \a other. + * + * \returns a reference to \c *this + * + * Example: \include MatrixBase_applyOnTheRight.cpp + * Output: \verbinclude MatrixBase_applyOnTheRight.out + */ +template +template +inline Derived& +MatrixBase::operator*=(const EigenBase &other) +{ + other.derived().applyThisOnTheRight(derived()); + return derived(); +} + +/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=(). + * + * Example: \include MatrixBase_applyOnTheRight.cpp + * Output: \verbinclude MatrixBase_applyOnTheRight.out + */ +template +template +inline void MatrixBase::applyOnTheRight(const EigenBase &other) +{ + other.derived().applyThisOnTheRight(derived()); +} + +/** replaces \c *this by \a other * \c *this. + * + * Example: \include MatrixBase_applyOnTheLeft.cpp + * Output: \verbinclude MatrixBase_applyOnTheLeft.out + */ +template +template +inline void MatrixBase::applyOnTheLeft(const EigenBase &other) +{ + other.derived().applyThisOnTheLeft(derived()); +} + } // end namespace Eigen #endif // EIGEN_MATRIXBASE_H diff --git a/otherlibs/eigen3/Eigen/src/Core/PermutationMatrix.h b/otherlibs/eigen3/Eigen/src/Core/PermutationMatrix.h index 4fc5dd3189..1297b8413f 100644 --- a/otherlibs/eigen3/Eigen/src/Core/PermutationMatrix.h +++ b/otherlibs/eigen3/Eigen/src/Core/PermutationMatrix.h @@ -553,7 +553,8 @@ struct permut_matrix_product_retval template inline void evalTo(Dest& dst) const { const Index n = Side==OnTheLeft ? rows() : cols(); - + // FIXME we need an is_same for expression that is not sensitive to constness. For instance + // is_same_xpr, Block >::value should be true. if(is_same::value && extract_data(dst) == extract_data(m_matrix)) { // apply the permutation inplace diff --git a/otherlibs/eigen3/Eigen/src/Core/PlainObjectBase.h b/otherlibs/eigen3/Eigen/src/Core/PlainObjectBase.h index af0a479c76..dd34b59e54 100644 --- a/otherlibs/eigen3/Eigen/src/Core/PlainObjectBase.h +++ b/otherlibs/eigen3/Eigen/src/Core/PlainObjectBase.h @@ -47,7 +47,10 @@ template<> struct check_rows_cols_for_overflow { } }; -template struct conservative_resize_like_impl; +template +struct conservative_resize_like_impl; template struct matrix_swap_impl; @@ -668,8 +671,10 @@ class PlainObjectBase : public internal::dense_xpr_base::type enum { ThisConstantIsPrivateInPlainObjectBase }; }; +namespace internal { + template -struct internal::conservative_resize_like_impl +struct conservative_resize_like_impl { typedef typename Derived::Index Index; static void run(DenseBase& _this, Index rows, Index cols) @@ -729,11 +734,14 @@ struct internal::conservative_resize_like_impl } }; -namespace internal { - +// Here, the specialization for vectors inherits from the general matrix case +// to allow calling .conservativeResize(rows,cols) on vectors. template struct conservative_resize_like_impl + : conservative_resize_like_impl { + using conservative_resize_like_impl::run; + typedef typename Derived::Index Index; static void run(DenseBase& _this, Index size) { diff --git a/otherlibs/eigen3/Eigen/src/Core/Ref.h b/otherlibs/eigen3/Eigen/src/Core/Ref.h index aba795bdb7..00d9e6d2b4 100644 --- a/otherlibs/eigen3/Eigen/src/Core/Ref.h +++ b/otherlibs/eigen3/Eigen/src/Core/Ref.h @@ -94,7 +94,8 @@ struct traits > typedef _PlainObjectType PlainObjectType; typedef _StrideType StrideType; enum { - Options = _Options + Options = _Options, + Flags = traits >::Flags | NestByRefBit }; template struct match { @@ -111,7 +112,7 @@ struct traits > }; typedef typename internal::conditional::type type; }; - + }; template diff --git a/otherlibs/eigen3/Eigen/src/Core/StableNorm.h b/otherlibs/eigen3/Eigen/src/Core/StableNorm.h index c83e955ee0..389d942753 100644 --- a/otherlibs/eigen3/Eigen/src/Core/StableNorm.h +++ b/otherlibs/eigen3/Eigen/src/Core/StableNorm.h @@ -17,16 +17,29 @@ namespace internal { template inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale) { - Scalar max = bl.cwiseAbs().maxCoeff(); - if (max>scale) + using std::max; + Scalar maxCoeff = bl.cwiseAbs().maxCoeff(); + + if (maxCoeff>scale) { - ssq = ssq * numext::abs2(scale/max); - scale = max; - invScale = Scalar(1)/scale; + ssq = ssq * numext::abs2(scale/maxCoeff); + Scalar tmp = Scalar(1)/maxCoeff; + if(tmp > NumTraits::highest()) + { + invScale = NumTraits::highest(); + scale = Scalar(1)/invScale; + } + else + { + scale = maxCoeff; + invScale = tmp; + } } - // TODO if the max is much much smaller than the current scale, + + // TODO if the maxCoeff is much much smaller than the current scale, // then we can neglect this sub vector - ssq += (bl*invScale).squaredNorm(); + if(scale>Scalar(0)) // if scale==0, then bl is 0 + ssq += (bl*invScale).squaredNorm(); } template diff --git a/otherlibs/eigen3/Eigen/src/Core/Transpose.h b/otherlibs/eigen3/Eigen/src/Core/Transpose.h index f21b3aa65f..22096ea2fd 100644 --- a/otherlibs/eigen3/Eigen/src/Core/Transpose.h +++ b/otherlibs/eigen3/Eigen/src/Core/Transpose.h @@ -284,7 +284,8 @@ struct inplace_transpose_selector { // non square matrix * Notice however that this method is only useful if you want to replace a matrix by its own transpose. * If you just need the transpose of a matrix, use transpose(). * - * \note if the matrix is not square, then \c *this must be a resizable matrix. + * \note if the matrix is not square, then \c *this must be a resizable matrix. + * This excludes (non-square) fixed-size matrices, block-expressions and maps. * * \sa transpose(), adjoint(), adjointInPlace() */ template @@ -315,6 +316,7 @@ inline void DenseBase::transposeInPlace() * If you just need the adjoint of a matrix, use adjoint(). * * \note if the matrix is not square, then \c *this must be a resizable matrix. + * This excludes (non-square) fixed-size matrices, block-expressions and maps. * * \sa transpose(), adjoint(), transposeInPlace() */ template diff --git a/otherlibs/eigen3/Eigen/src/Core/VectorwiseOp.h b/otherlibs/eigen3/Eigen/src/Core/VectorwiseOp.h index 511564875a..d5ab036647 100644 --- a/otherlibs/eigen3/Eigen/src/Core/VectorwiseOp.h +++ b/otherlibs/eigen3/Eigen/src/Core/VectorwiseOp.h @@ -50,7 +50,7 @@ struct traits > MaxColsAtCompileTime = Direction==Horizontal ? 1 : MatrixType::MaxColsAtCompileTime, Flags0 = (unsigned int)_MatrixTypeNested::Flags & HereditaryBits, Flags = (Flags0 & ~RowMajorBit) | (RowsAtCompileTime == 1 ? RowMajorBit : 0), - TraversalSize = Direction==Vertical ? RowsAtCompileTime : ColsAtCompileTime + TraversalSize = Direction==Vertical ? MatrixType::RowsAtCompileTime : MatrixType::ColsAtCompileTime }; #if EIGEN_GNUC_AT_LEAST(3,4) typedef typename MemberOp::template Cost CostOpType; @@ -58,7 +58,8 @@ struct traits > typedef typename MemberOp::template Cost CostOpType; #endif enum { - CoeffReadCost = TraversalSize * traits<_MatrixTypeNested>::CoeffReadCost + int(CostOpType::value) + CoeffReadCost = TraversalSize==Dynamic ? Dynamic + : TraversalSize * traits<_MatrixTypeNested>::CoeffReadCost + int(CostOpType::value) }; }; } diff --git a/otherlibs/eigen3/Eigen/src/Core/arch/SSE/MathFunctions.h b/otherlibs/eigen3/Eigen/src/Core/arch/SSE/MathFunctions.h index 3376a984e5..99cbd0d95b 100644 --- a/otherlibs/eigen3/Eigen/src/Core/arch/SSE/MathFunctions.h +++ b/otherlibs/eigen3/Eigen/src/Core/arch/SSE/MathFunctions.h @@ -442,8 +442,11 @@ Packet4f pcos(const Packet4f& _x) return _mm_xor_ps(y, sign_bit); } +#if EIGEN_FAST_MATH + // This is based on Quake3's fast inverse square root. // For detail see here: http://www.beyond3d.com/content/articles/8/ +// It lacks 1 (or 2 bits in some rare cases) of precision, and does not handle negative, +inf, or denormalized numbers correctly. template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f psqrt(const Packet4f& _x) { @@ -457,6 +460,14 @@ Packet4f psqrt(const Packet4f& _x) return pmul(_x,x); } +#else + +template<> EIGEN_STRONG_INLINE Packet4f psqrt(const Packet4f& x) { return _mm_sqrt_ps(x); } + +#endif + +template<> EIGEN_STRONG_INLINE Packet2d psqrt(const Packet2d& x) { return _mm_sqrt_pd(x); } + } // end namespace internal } // end namespace Eigen diff --git a/otherlibs/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h b/otherlibs/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h index e256f4bac4..fc8ae50fed 100644 --- a/otherlibs/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/otherlibs/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h @@ -83,7 +83,8 @@ template<> struct packet_traits : default_packet_traits size=2, HasDiv = 1, - HasExp = 1 + HasExp = 1, + HasSqrt = 1 }; }; template<> struct packet_traits : default_packet_traits @@ -507,8 +508,8 @@ template<> EIGEN_STRONG_INLINE int predux_min(const Packet4i& a) // for GCC (eg., it does not like using std::min after the pstore !!) EIGEN_ALIGN16 int aux[4]; pstore(aux, a); - register int aux0 = aux[0] EIGEN_STRONG_INLINE int predux_max(const Packet4i& a) // for GCC (eg., it does not like using std::min after the pstore !!) EIGEN_ALIGN16 int aux[4]; pstore(aux, a); - register int aux0 = aux[0]>aux[1] ? aux[0] : aux[1]; - register int aux2 = aux[2]>aux[3] ? aux[2] : aux[3]; + int aux0 = aux[0]>aux[1] ? aux[0] : aux[1]; + int aux2 = aux[2]>aux[3] ? aux[2] : aux[3]; return aux0>aux2 ? aux0 : aux2; } diff --git a/otherlibs/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/otherlibs/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 780fa74d3f..bcdca5b0d2 100644 --- a/otherlibs/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/otherlibs/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -1128,6 +1128,8 @@ EIGEN_DONT_INLINE void gemm_pack_lhs::size }; EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK LHS"); + EIGEN_UNUSED_VARIABLE(stride) + EIGEN_UNUSED_VARIABLE(offset) eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride)); eigen_assert( (StorageOrder==RowMajor) || ((Pack1%PacketSize)==0 && Pack1<=4*PacketSize) ); conj_if::IsComplex && Conjugate> cj; @@ -1215,6 +1217,8 @@ EIGEN_DONT_INLINE void gemm_pack_rhs=depth && offset<=stride)); conj_if::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; @@ -1266,6 +1270,8 @@ EIGEN_DONT_INLINE void gemm_pack_rhs=depth && offset<=stride)); conj_if::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; diff --git a/otherlibs/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h b/otherlibs/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h index c1cb784988..09387703ef 100644 --- a/otherlibs/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/otherlibs/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h @@ -52,11 +52,7 @@ EIGEN_DONT_INLINE static void run( Index rows, Index cols, const LhsScalar* lhs, Index lhsStride, const RhsScalar* rhs, Index rhsIncr, - ResScalar* res, Index - #ifdef EIGEN_INTERNAL_DEBUGGING - resIncr - #endif - , RhsScalar alpha); + ResScalar* res, Index resIncr, RhsScalar alpha); }; template @@ -64,12 +60,9 @@ EIGEN_DONT_INLINE void general_matrix_vector_product(&lhs0[i]), ptmp0, pload(&res[i]))); + pstore(&res[i], pcj.pmadd(pload(&lhs0[i]), ptmp0, pload(&res[i]))); else for (Index i = alignedStart;i(&lhs0[i]), ptmp0, pload(&res[i]))); diff --git a/otherlibs/eigen3/Eigen/src/Core/products/SelfadjointMatrixVector.h b/otherlibs/eigen3/Eigen/src/Core/products/SelfadjointMatrixVector.h index c40e80f53c..f698f67f9b 100644 --- a/otherlibs/eigen3/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/otherlibs/eigen3/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -79,8 +79,8 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product(t0); @@ -147,7 +147,7 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_productx || (EIGEN_WORLD_VERSION>=x && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ @@ -238,7 +238,12 @@ #endif // Suppresses 'unused variable' warnings. -#define EIGEN_UNUSED_VARIABLE(var) (void)var; +namespace Eigen { + namespace internal { + template void ignore_unused_variable(const T&) {} + } +} +#define EIGEN_UNUSED_VARIABLE(var) Eigen::internal::ignore_unused_variable(var); #if !defined(EIGEN_ASM_COMMENT) #if (defined __GNUC__) && ( defined(__i386__) || defined(__x86_64__) ) diff --git a/otherlibs/eigen3/Eigen/src/Core/util/Memory.h b/otherlibs/eigen3/Eigen/src/Core/util/Memory.h index 451535a0c8..cacbd02fd1 100644 --- a/otherlibs/eigen3/Eigen/src/Core/util/Memory.h +++ b/otherlibs/eigen3/Eigen/src/Core/util/Memory.h @@ -578,7 +578,7 @@ template class aligned_stack_memory_handler */ #ifdef EIGEN_ALLOCA - #ifdef __arm__ + #if defined(__arm__) || defined(_WIN32) #define EIGEN_ALIGNED_ALLOCA(SIZE) reinterpret_cast((reinterpret_cast(EIGEN_ALLOCA(SIZE+16)) & ~(size_t(15))) + 16) #else #define EIGEN_ALIGNED_ALLOCA EIGEN_ALLOCA @@ -634,7 +634,9 @@ template class aligned_stack_memory_handler /* memory allocated we can safely let the default implementation handle */ \ /* this particular case. */ \ static void *operator new(size_t size, void *ptr) { return ::operator new(size,ptr); } \ + static void *operator new[](size_t size, void* ptr) { return ::operator new[](size,ptr); } \ void operator delete(void * memory, void *ptr) throw() { return ::operator delete(memory,ptr); } \ + void operator delete[](void * memory, void *ptr) throw() { return ::operator delete[](memory,ptr); } \ /* nothrow-new (returns zero instead of std::bad_alloc) */ \ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \ void operator delete(void *ptr, const std::nothrow_t&) throw() { \ @@ -729,15 +731,6 @@ class aligned_allocator ::new( p ) T( value ); } - // Support for c++11 -#if (__cplusplus >= 201103L) - template - void construct(pointer p, Args&&... args) - { - ::new(p) T(std::forward(args)...); - } -#endif - void destroy( pointer p ) { p->~T(); diff --git a/otherlibs/eigen3/Eigen/src/Eigen2Support/SVD.h b/otherlibs/eigen3/Eigen/src/Eigen2Support/SVD.h index 077d26d543..3d03d2288d 100644 --- a/otherlibs/eigen3/Eigen/src/Eigen2Support/SVD.h +++ b/otherlibs/eigen3/Eigen/src/Eigen2Support/SVD.h @@ -512,8 +512,7 @@ template template bool SVD::solve(const MatrixBase &b, ResultType* result) const { - const int rows = m_matU.rows(); - ei_assert(b.rows() == rows); + ei_assert(b.rows() == m_matU.rows()); Scalar maxVal = m_sigma.cwise().abs().maxCoeff(); for (int j=0; j::squaredExteriorDistance(const Align return dist2; } -/** \defgroup alignedboxtypedefs Global aligned box typedefs - * \ingroup eigen_grp +/** \defgroup alignedboxtypedefs Global aligned box typedefs * * \ingroup Geometry_Module * diff --git a/otherlibs/eigen3/Eigen/src/Geometry/EulerAngles.h b/otherlibs/eigen3/Eigen/src/Geometry/EulerAngles.h index 97984d590e..82802fb43c 100644 --- a/otherlibs/eigen3/Eigen/src/Geometry/EulerAngles.h +++ b/otherlibs/eigen3/Eigen/src/Geometry/EulerAngles.h @@ -28,7 +28,7 @@ namespace Eigen { * * AngleAxisf(ea[2], Vector3f::UnitZ()); \endcode * This corresponds to the right-multiply conventions (with right hand side frames). * - * The returned angles are in the ranges [0:pi]x[0:pi]x[-pi:pi]. + * The returned angles are in the ranges [0:pi]x[-pi:pi]x[-pi:pi]. * * \sa class AngleAxis */ diff --git a/otherlibs/eigen3/Eigen/src/Geometry/Quaternion.h b/otherlibs/eigen3/Eigen/src/Geometry/Quaternion.h index e135f2b663..9fee0c9198 100644 --- a/otherlibs/eigen3/Eigen/src/Geometry/Quaternion.h +++ b/otherlibs/eigen3/Eigen/src/Geometry/Quaternion.h @@ -150,10 +150,6 @@ class QuaternionBase : public RotationBase /** \returns the conjugated quaternion */ Quaternion conjugate() const; - /** \returns an interpolation for a constant motion between \a other and \c *this - * \a t in [0;1] - * see http://en.wikipedia.org/wiki/Slerp - */ template Quaternion slerp(const Scalar& t, const QuaternionBase& other) const; /** \returns \c true if \c *this is approximately equal to \a other, within the precision @@ -194,11 +190,11 @@ class QuaternionBase : public RotationBase * \brief The quaternion class used to represent 3D orientations and rotations * * \tparam _Scalar the scalar type, i.e., the type of the coefficients - * \tparam _Options controls the memory alignement of the coeffecients. Can be \# AutoAlign or \# DontAlign. Default is AutoAlign. + * \tparam _Options controls the memory alignment of the coefficients. Can be \# AutoAlign or \# DontAlign. Default is AutoAlign. * * This class represents a quaternion \f$ w+xi+yj+zk \f$ that is a convenient representation of * orientations and rotations of objects in three dimensions. Compared to other representations - * like Euler angles or 3x3 matrices, quatertions offer the following advantages: + * like Euler angles or 3x3 matrices, quaternions offer the following advantages: * \li \b compact storage (4 scalars) * \li \b efficient to compose (28 flops), * \li \b stable spherical interpolation @@ -385,7 +381,7 @@ class Map, _Options > /** Constructs a Mapped Quaternion object from the pointer \a coeffs * - * The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order: + * The pointer \a coeffs must reference the four coefficients of Quaternion in the following order: * \code *coeffs == {x, y, z, w} \endcode * * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */ @@ -399,16 +395,16 @@ class Map, _Options > }; /** \ingroup Geometry_Module - * Map an unaligned array of single precision scalar as a quaternion */ + * Map an unaligned array of single precision scalars as a quaternion */ typedef Map, 0> QuaternionMapf; /** \ingroup Geometry_Module - * Map an unaligned array of double precision scalar as a quaternion */ + * Map an unaligned array of double precision scalars as a quaternion */ typedef Map, 0> QuaternionMapd; /** \ingroup Geometry_Module - * Map a 16-bits aligned array of double precision scalars as a quaternion */ + * Map a 16-byte aligned array of single precision scalars as a quaternion */ typedef Map, Aligned> QuaternionMapAlignedf; /** \ingroup Geometry_Module - * Map a 16-bits aligned array of double precision scalars as a quaternion */ + * Map a 16-byte aligned array of double precision scalars as a quaternion */ typedef Map, Aligned> QuaternionMapAlignedd; /*************************************************************************** @@ -579,7 +575,7 @@ inline Derived& QuaternionBase::setFromTwoVectors(const MatrixBase accuraletly compute the rotation axis by computing the + // => accurately compute the rotation axis by computing the // intersection of the two planes. This is done by solving: // x^T v0 = 0 // x^T v1 = 0 @@ -677,8 +673,13 @@ QuaternionBase::angularDistance(const QuaternionBase& oth return static_cast(2 * acos(d)); } + + /** \returns the spherical linear interpolation between the two quaternions - * \c *this and \a other at the parameter \a t + * \c *this and \a other at the parameter \a t in [0;1]. + * + * This represents an interpolation for a constant motion between \c *this and \a other, + * see also http://en.wikipedia.org/wiki/Slerp. */ template template diff --git a/otherlibs/eigen3/Eigen/src/Geometry/Transform.h b/otherlibs/eigen3/Eigen/src/Geometry/Transform.h index 887e718d67..498fea41a9 100644 --- a/otherlibs/eigen3/Eigen/src/Geometry/Transform.h +++ b/otherlibs/eigen3/Eigen/src/Geometry/Transform.h @@ -530,9 +530,9 @@ class Transform inline Transform& operator=(const UniformScaling& t); inline Transform& operator*=(const UniformScaling& s) { return scale(s.factor()); } - inline Transform operator*(const UniformScaling& s) const + inline Transform operator*(const UniformScaling& s) const { - Transform res = *this; + Transform res = *this; res.scale(s.factor()); return res; } diff --git a/otherlibs/eigen3/Eigen/src/QR/ColPivHouseholderQR.h b/otherlibs/eigen3/Eigen/src/QR/ColPivHouseholderQR.h index 8b01f8179e..bec85810cc 100644 --- a/otherlibs/eigen3/Eigen/src/QR/ColPivHouseholderQR.h +++ b/otherlibs/eigen3/Eigen/src/QR/ColPivHouseholderQR.h @@ -349,7 +349,7 @@ template class ColPivHouseholderQR return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. - : NumTraits::epsilon() * m_qr.diagonalSize(); + : NumTraits::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. diff --git a/otherlibs/eigen3/Eigen/src/QR/FullPivHouseholderQR.h b/otherlibs/eigen3/Eigen/src/QR/FullPivHouseholderQR.h index 0dd5ad3472..6168e7abfb 100644 --- a/otherlibs/eigen3/Eigen/src/QR/FullPivHouseholderQR.h +++ b/otherlibs/eigen3/Eigen/src/QR/FullPivHouseholderQR.h @@ -63,9 +63,10 @@ template class FullPivHouseholderQR typedef typename MatrixType::Index Index; typedef internal::FullPivHouseholderQRMatrixQReturnType MatrixQReturnType; typedef typename internal::plain_diag_type::type HCoeffsType; - typedef Matrix IntRowVectorType; + typedef Matrix IntDiagSizeVectorType; typedef PermutationMatrix PermutationType; - typedef typename internal::plain_col_type::type IntColVectorType; typedef typename internal::plain_row_type::type RowVectorType; typedef typename internal::plain_col_type::type ColVectorType; @@ -93,10 +94,10 @@ template class FullPivHouseholderQR FullPivHouseholderQR(Index rows, Index cols) : m_qr(rows, cols), m_hCoeffs((std::min)(rows,cols)), - m_rows_transpositions(rows), - m_cols_transpositions(cols), + m_rows_transpositions((std::min)(rows,cols)), + m_cols_transpositions((std::min)(rows,cols)), m_cols_permutation(cols), - m_temp((std::min)(rows,cols)), + m_temp(cols), m_isInitialized(false), m_usePrescribedThreshold(false) {} @@ -115,10 +116,10 @@ template class FullPivHouseholderQR FullPivHouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), - m_rows_transpositions(matrix.rows()), - m_cols_transpositions(matrix.cols()), + m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), + m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), m_cols_permutation(matrix.cols()), - m_temp((std::min)(matrix.rows(), matrix.cols())), + m_temp(matrix.cols()), m_isInitialized(false), m_usePrescribedThreshold(false) { @@ -126,11 +127,12 @@ template class FullPivHouseholderQR } /** This method finds a solution x to the equation Ax=b, where A is the matrix of which - * *this is the QR decomposition, if any exists. + * \c *this is the QR decomposition. * * \param b the right-hand-side of the equation to solve. * - * \returns a solution. + * \returns the exact or least-square solution if the rank is greater or equal to the number of columns of A, + * and an arbitrary solution otherwise. * * \note The case where b is a matrix is not yet implemented. Also, this * code is space inefficient. @@ -172,7 +174,7 @@ template class FullPivHouseholderQR } /** \returns a const reference to the vector of indices representing the rows transpositions */ - const IntColVectorType& rowsTranspositions() const + const IntDiagSizeVectorType& rowsTranspositions() const { eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); return m_rows_transpositions; @@ -344,7 +346,7 @@ template class FullPivHouseholderQR return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. - : NumTraits::epsilon() * m_qr.diagonalSize(); + : NumTraits::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. @@ -368,8 +370,8 @@ template class FullPivHouseholderQR protected: MatrixType m_qr; HCoeffsType m_hCoeffs; - IntColVectorType m_rows_transpositions; - IntRowVectorType m_cols_transpositions; + IntDiagSizeVectorType m_rows_transpositions; + IntDiagSizeVectorType m_cols_transpositions; PermutationType m_cols_permutation; RowVectorType m_temp; bool m_isInitialized, m_usePrescribedThreshold; @@ -415,10 +417,10 @@ FullPivHouseholderQR& FullPivHouseholderQR::compute(cons m_temp.resize(cols); - m_precision = NumTraits::epsilon() * size; + m_precision = NumTraits::epsilon() * RealScalar(size); - m_rows_transpositions.resize(matrix.rows()); - m_cols_transpositions.resize(matrix.cols()); + m_rows_transpositions.resize(size); + m_cols_transpositions.resize(size); Index number_of_transpositions = 0; RealScalar biggest(0); @@ -516,17 +518,6 @@ struct solve_retval, Rhs> dec().hCoeffs().coeff(k), &temp.coeffRef(0)); } - if(!dec().isSurjective()) - { - // is c is in the image of R ? - RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff(); - RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff(); - // FIXME brain dead - const RealScalar m_precision = NumTraits::epsilon() * (std::min)(rows,cols); - // this internal:: prefix is needed by at least gcc 3.4 and ICC - if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision)) - return; - } dec().matrixQR() .topLeftCorner(dec().rank(), dec().rank()) .template triangularView() @@ -548,14 +539,14 @@ template struct FullPivHouseholderQRMatrixQReturnType { public: typedef typename MatrixType::Index Index; - typedef typename internal::plain_col_type::type IntColVectorType; + typedef typename FullPivHouseholderQR::IntDiagSizeVectorType IntDiagSizeVectorType; typedef typename internal::plain_diag_type::type HCoeffsType; typedef Matrix WorkVectorType; FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr, const HCoeffsType& hCoeffs, - const IntColVectorType& rowsTranspositions) + const IntDiagSizeVectorType& rowsTranspositions) : m_qr(qr), m_hCoeffs(hCoeffs), m_rowsTranspositions(rowsTranspositions) @@ -595,7 +586,7 @@ template struct FullPivHouseholderQRMatrixQReturnType protected: typename MatrixType::Nested m_qr; typename HCoeffsType::Nested m_hCoeffs; - typename IntColVectorType::Nested m_rowsTranspositions; + typename IntDiagSizeVectorType::Nested m_rowsTranspositions; }; } // end namespace internal diff --git a/otherlibs/eigen3/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/otherlibs/eigen3/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h index aa41f434c0..a2cc2a9e26 100644 --- a/otherlibs/eigen3/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h +++ b/otherlibs/eigen3/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h @@ -64,7 +64,8 @@ class SPQR typedef PermutationMatrix PermutationType; public: SPQR() - : m_ordering(SPQR_ORDERING_DEFAULT), + : m_isInitialized(false), + m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits::epsilon()) { @@ -72,7 +73,8 @@ class SPQR } SPQR(const _MatrixType& matrix) - : m_ordering(SPQR_ORDERING_DEFAULT), + : m_isInitialized(false), + m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits::epsilon()) { @@ -82,16 +84,22 @@ class SPQR ~SPQR() { - // Calls SuiteSparseQR_free() + SPQR_free(); + cholmod_l_finish(&m_cc); + } + void SPQR_free() + { cholmod_l_free_sparse(&m_H, &m_cc); cholmod_l_free_sparse(&m_cR, &m_cc); cholmod_l_free_dense(&m_HTau, &m_cc); std::free(m_E); std::free(m_HPinv); - cholmod_l_finish(&m_cc); } + void compute(const _MatrixType& matrix) { + if(m_isInitialized) SPQR_free(); + MatrixType mat(matrix); cholmod_sparse A; A = viewAsCholmod(mat); @@ -139,7 +147,7 @@ class SPQR eigen_assert(b.cols()==1 && "This method is for vectors only"); //Compute Q^T * b - Dest y; + typename Dest::PlainObject y; y = matrixQ().transpose() * b; // Solves with the triangular matrix R Index rk = this->rank(); diff --git a/otherlibs/eigen3/Eigen/src/SVD/JacobiSVD.h b/otherlibs/eigen3/Eigen/src/SVD/JacobiSVD.h index e77396adf2..f44995cd39 100644 --- a/otherlibs/eigen3/Eigen/src/SVD/JacobiSVD.h +++ b/otherlibs/eigen3/Eigen/src/SVD/JacobiSVD.h @@ -10,7 +10,7 @@ #ifndef EIGEN_JACOBISVD_H #define EIGEN_JACOBISVD_H -namespace Eigen { +namespace Eigen { namespace internal { // forward declaration (needed by ICC) @@ -380,7 +380,10 @@ struct svd_precondition_2x2_block_to_be_real z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.row(p) *= z; if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); - z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + if(work_matrix.coeff(q,q)!=Scalar(0)) + z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + else + z = Scalar(0); work_matrix.row(q) *= z; if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); } @@ -666,19 +669,6 @@ template class JacobiSVD inline Index rows() const { return m_rows; } inline Index cols() const { return m_cols; } - void pinv( MatrixType& pinvmat) const // New method added from outside - { - eigen_assert(m_isInitialized && "SVD is not initialized."); - double pinvtoler=1.e-6; // choose your tolerance wisely! - SingularValuesType singularValues_inv=m_singularValues; - for ( long i=0; i pinvtoler ) - singularValues_inv(i)=1.0/m_singularValues(i); - else singularValues_inv(i)=0; - } - pinvmat= (m_matrixV*singularValues_inv.asDiagonal()*m_matrixU.transpose()); - } - private: void allocate(Index rows, Index cols, unsigned int computationOptions); @@ -745,7 +735,7 @@ void JacobiSVD::allocate(Index rows, Index cols, u : m_computeThinV ? m_diagSize : 0); m_workMatrix.resize(m_diagSize, m_diagSize); - + if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); } @@ -865,7 +855,7 @@ struct solve_retval, Rhs> Matrix tmp; Index nonzeroSingVals = dec().nonzeroSingularValues(); - + tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs(); tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp; dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp; diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/otherlibs/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h index 4b13f08d4e..5c320e2d2d 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h @@ -66,9 +66,9 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r } // unordered insertion - for(int k=0; k1) std::sort(indices.data(),indices.data()+nnz); - for(int k=0; k RowMajorMatrix; - typedef SparseMatrix ColMajorMatrix; + typedef SparseMatrix RowMajorMatrix; + typedef SparseMatrix ColMajorMatrix; ColMajorMatrix resCol(lhs.rows(),rhs.cols()); internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); // sort the non zeros: @@ -149,7 +149,7 @@ struct conservative_sparse_sparse_product_selector RowMajorMatrix; + typedef SparseMatrix RowMajorMatrix; RowMajorMatrix rhsRow = rhs; RowMajorMatrix resRow(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl(rhsRow, lhs, resRow); @@ -162,7 +162,7 @@ struct conservative_sparse_sparse_product_selector RowMajorMatrix; + typedef SparseMatrix RowMajorMatrix; RowMajorMatrix lhsRow = lhs; RowMajorMatrix resRow(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl(rhs, lhsRow, resRow); @@ -175,7 +175,7 @@ struct conservative_sparse_sparse_product_selector RowMajorMatrix; + typedef SparseMatrix RowMajorMatrix; RowMajorMatrix resRow(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl(rhs, lhs, resRow); res = resRow; @@ -190,7 +190,7 @@ struct conservative_sparse_sparse_product_selector ColMajorMatrix; + typedef SparseMatrix ColMajorMatrix; ColMajorMatrix resCol(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); res = resCol; @@ -202,7 +202,7 @@ struct conservative_sparse_sparse_product_selector ColMajorMatrix; + typedef SparseMatrix ColMajorMatrix; ColMajorMatrix lhsCol = lhs; ColMajorMatrix resCol(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl(lhsCol, rhs, resCol); @@ -215,7 +215,7 @@ struct conservative_sparse_sparse_product_selector ColMajorMatrix; + typedef SparseMatrix ColMajorMatrix; ColMajorMatrix rhsCol = rhs; ColMajorMatrix resCol(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl(lhs, rhsCol, resCol); @@ -228,8 +228,8 @@ struct conservative_sparse_sparse_product_selector RowMajorMatrix; - typedef SparseMatrix ColMajorMatrix; + typedef SparseMatrix RowMajorMatrix; + typedef SparseMatrix ColMajorMatrix; RowMajorMatrix resRow(lhs.rows(),rhs.cols()); internal::conservative_sparse_sparse_product_impl(rhs, lhs, resRow); // sort the non zeros: diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/MappedSparseMatrix.h b/otherlibs/eigen3/Eigen/src/SparseCore/MappedSparseMatrix.h index 93cd4832de..ab1a266a90 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/MappedSparseMatrix.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/MappedSparseMatrix.h @@ -50,6 +50,8 @@ class MappedSparseMatrix inline Index cols() const { return IsRowMajor ? m_innerSize : m_outerSize; } inline Index innerSize() const { return m_innerSize; } inline Index outerSize() const { return m_outerSize; } + + bool isCompressed() const { return true; } //---------------------------------------- // direct access interface diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/SparseBlock.h b/otherlibs/eigen3/Eigen/src/SparseCore/SparseBlock.h index 0b3e193dbc..16a20a5744 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/SparseBlock.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/SparseBlock.h @@ -66,6 +66,8 @@ class BlockImpl typename XprType::Nested m_matrix; Index m_outerStart; const internal::variable_if_dynamic m_outerSize; + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) }; @@ -391,6 +393,8 @@ class BlockImpl protected: friend class InnerIterator; friend class ReverseInnerIterator; + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) typename XprType::Nested m_matrix; const internal::variable_if_dynamic m_startRow; diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/SparseColEtree.h b/otherlibs/eigen3/Eigen/src/SparseCore/SparseColEtree.h index f89ca38144..f8745f4610 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/SparseColEtree.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/SparseColEtree.h @@ -63,6 +63,7 @@ int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowEl typedef typename MatrixType::Index Index; Index nc = mat.cols(); // Number of columns Index m = mat.rows(); + Index diagSize = (std::min)(nc,m); IndexVector root(nc); // root of subtree of etree root.setZero(); IndexVector pp(nc); // disjoint sets @@ -72,7 +73,7 @@ int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowEl Index row,col; firstRowElt.resize(m); firstRowElt.setConstant(nc); - firstRowElt.segment(0, nc).setLinSpaced(nc, 0, nc-1); + firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1); bool found_diag; for (col = 0; col < nc; col++) { @@ -91,7 +92,7 @@ int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowEl Index rset, cset, rroot; for (col = 0; col < nc; col++) { - found_diag = false; + found_diag = col>=m; pp(col) = col; cset = col; root(cset) = col; @@ -105,6 +106,7 @@ int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowEl Index i = col; if(it) i = it.index(); if (i == col) found_diag = true; + row = firstRowElt(i); if (row >= col) continue; rset = internal::etree_find(row, pp); // Find the name of the set containing row diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/SparseDenseProduct.h b/otherlibs/eigen3/Eigen/src/SparseCore/SparseDenseProduct.h index 30975c29cc..54fd633a10 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/SparseDenseProduct.h @@ -125,7 +125,7 @@ class SparseDenseOuterProduct::InnerIterator : public _LhsNes inline Scalar value() const { return Base::value() * m_factor; } protected: - int m_outer; + Index m_outer; Scalar m_factor; }; @@ -155,7 +155,7 @@ struct sparse_time_dense_product_impl::Constant(outerSize(), 2)); } return insertUncompressed(row,col); } @@ -402,7 +402,7 @@ class SparseMatrix * \sa insertBack, insertBackByOuterInner */ inline void startVec(Index outer) { - eigen_assert(m_outerIndex[outer]==int(m_data.size()) && "You must call startVec for each inner vector sequentially"); + eigen_assert(m_outerIndex[outer]==Index(m_data.size()) && "You must call startVec for each inner vector sequentially"); eigen_assert(m_outerIndex[outer+1]==0 && "You must call startVec for each inner vector sequentially"); m_outerIndex[outer+1] = m_outerIndex[outer]; } @@ -480,7 +480,7 @@ class SparseMatrix if(m_innerNonZeros != 0) return; m_innerNonZeros = static_cast(std::malloc(m_outerSize * sizeof(Index))); - for (int i = 0; i < m_outerSize; i++) + for (Index i = 0; i < m_outerSize; i++) { m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; } @@ -752,8 +752,8 @@ class SparseMatrix else for (Index i=0; i trMat(mat.rows(),mat.cols()); - if(begin wi(trMat.outerSize()); wi.setZero(); for(InputIterator it(begin); it!=end; ++it) { @@ -1018,11 +1019,11 @@ void SparseMatrix::sumupDuplicates() { eigen_assert(!isCompressed()); // TODO, in practice we should be able to use m_innerNonZeros for that task - VectorXi wi(innerSize()); + Matrix wi(innerSize()); wi.fill(-1); Index count = 0; // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers - for(int j=0; j& SparseMatrix positions(dest.outerSize()); for (Index j=0; j class SparseMatrixBase : public EigenBase } else { - SparseMatrix trans = m; - s << static_cast >&>(trans); + SparseMatrix trans = m; + s << static_cast >&>(trans); } } return s; diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/SparsePermutation.h b/otherlibs/eigen3/Eigen/src/SparseCore/SparsePermutation.h index b897b7595b..b85be93f6f 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/SparsePermutation.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/SparsePermutation.h @@ -57,7 +57,7 @@ struct permut_sparsematrix_product_retval if(MoveOuter) { SparseMatrix tmp(m_matrix.rows(), m_matrix.cols()); - VectorXi sizes(m_matrix.outerSize()); + Matrix sizes(m_matrix.outerSize()); for(Index j=0; j tmp(m_matrix.rows(), m_matrix.cols()); - VectorXi sizes(tmp.outerSize()); + Matrix sizes(tmp.outerSize()); sizes.setZero(); PermutationMatrix perm; if((Side==OnTheLeft) ^ Transposed) diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/SparseProduct.h b/otherlibs/eigen3/Eigen/src/SparseCore/SparseProduct.h index 70b6480efe..cf76630700 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/SparseProduct.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/SparseProduct.h @@ -16,6 +16,7 @@ template struct SparseSparseProductReturnType { typedef typename internal::traits::Scalar Scalar; + typedef typename internal::traits::Index Index; enum { LhsRowMajor = internal::traits::Flags & RowMajorBit, RhsRowMajor = internal::traits::Flags & RowMajorBit, @@ -24,11 +25,11 @@ struct SparseSparseProductReturnType }; typedef typename internal::conditional, + SparseMatrix, typename internal::nested::type>::type LhsNested; typedef typename internal::conditional, + SparseMatrix, typename internal::nested::type>::type RhsNested; typedef SparseSparseProduct Type; diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/SparseSparseProductWithPruning.h b/otherlibs/eigen3/Eigen/src/SparseCore/SparseSparseProductWithPruning.h index 70857c7b6e..fcc18f5c9c 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/SparseSparseProductWithPruning.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/SparseSparseProductWithPruning.h @@ -27,7 +27,7 @@ static void sparse_sparse_product_with_pruning_impl(const Lhs& lhs, const Rhs& r // make sure to call innerSize/outerSize since we fake the storage order. Index rows = lhs.innerSize(); Index cols = rhs.outerSize(); - //int size = lhs.outerSize(); + //Index size = lhs.outerSize(); eigen_assert(lhs.outerSize() == rhs.innerSize()); // allocate a temporary buffer @@ -100,7 +100,7 @@ struct sparse_sparse_product_with_pruning_selector SparseTemporaryType; + typedef SparseMatrix SparseTemporaryType; SparseTemporaryType _res(res.rows(), res.cols()); internal::sparse_sparse_product_with_pruning_impl(lhs, rhs, _res, tolerance); res = _res; @@ -126,10 +126,11 @@ struct sparse_sparse_product_with_pruning_selector ColMajorMatrix; - ColMajorMatrix colLhs(lhs); - ColMajorMatrix colRhs(rhs); - internal::sparse_sparse_product_with_pruning_impl(colLhs, colRhs, res, tolerance); + typedef SparseMatrix ColMajorMatrixLhs; + typedef SparseMatrix ColMajorMatrixRhs; + ColMajorMatrixLhs colLhs(lhs); + ColMajorMatrixRhs colRhs(rhs); + internal::sparse_sparse_product_with_pruning_impl(colLhs, colRhs, res, tolerance); // let's transpose the product to get a column x column product // typedef SparseMatrix SparseTemporaryType; diff --git a/otherlibs/eigen3/Eigen/src/SparseCore/SparseUtil.h b/otherlibs/eigen3/Eigen/src/SparseCore/SparseUtil.h index 064a407072..05023858b1 100644 --- a/otherlibs/eigen3/Eigen/src/SparseCore/SparseUtil.h +++ b/otherlibs/eigen3/Eigen/src/SparseCore/SparseUtil.h @@ -143,7 +143,7 @@ template struct plain_matrix_type * * \sa SparseMatrix::setFromTriplets() */ -template +template::Index > class Triplet { public: diff --git a/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU.h b/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU.h index dd9eab2c2c..1d592f2c8c 100644 --- a/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU.h +++ b/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU.h @@ -219,9 +219,9 @@ class SparseLU : public internal::SparseLUImpl - bool _solve(const MatrixBase &B, MatrixBase &_X) const + bool _solve(const MatrixBase &B, MatrixBase &X_base) const { - Dest& X(_X.derived()); + Dest& X(X_base.derived()); eigen_assert(m_factorizationIsOk && "The matrix should be factorized first"); EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); @@ -229,8 +229,10 @@ class SparseLU : public internal::SparseLUImplmatrixL().solveInPlace(X); diff --git a/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU_Memory.h b/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU_Memory.h index a5158025c5..1ffa7d54e9 100644 --- a/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/otherlibs/eigen3/Eigen/src/SparseLU/SparseLU_Memory.h @@ -70,23 +70,30 @@ Index SparseLUImpl::expand(VectorType& vec, Index& length, Index if(num_expansions == 0 || keep_prev) new_len = length ; // First time allocate requested else - new_len = Index(alpha * length); + new_len = (std::max)(length+1,Index(alpha * length)); VectorType old_vec; // Temporary vector to hold the previous values if (nbElts > 0 ) old_vec = vec.segment(0,nbElts); //Allocate or expand the current vector - try +#ifdef EIGEN_EXCEPTIONS + try +#endif { vec.resize(new_len); } +#ifdef EIGEN_EXCEPTIONS catch(std::bad_alloc& ) +#else + if(!vec.size()) +#endif { - if ( !num_expansions ) + if (!num_expansions) { // First time to allocate from LUMemInit() - throw; // Pass the exception to LUMemInit() which has a try... catch block + // Let LUMemInit() deals with it. + return -1; } if (keep_prev) { @@ -100,12 +107,18 @@ Index SparseLUImpl::expand(VectorType& vec, Index& length, Index do { alpha = (alpha + 1)/2; - new_len = Index(alpha * length); + new_len = (std::max)(length+1,Index(alpha * length)); +#ifdef EIGEN_EXCEPTIONS try +#endif { vec.resize(new_len); } +#ifdef EIGEN_EXCEPTIONS catch(std::bad_alloc& ) +#else + if (!vec.size()) +#endif { tries += 1; if ( tries > 10) return new_len; @@ -139,10 +152,9 @@ template Index SparseLUImpl::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t& glu) { Index& num_expansions = glu.num_expansions; //No memory expansions so far - num_expansions = 0; - glu.nzumax = glu.nzlumax = (std::max)(fillratio * annz, m*n); // estimated number of nonzeros in U + num_expansions = 0; + glu.nzumax = glu.nzlumax = (std::min)(fillratio * annz / n, m) * n; // estimated number of nonzeros in U glu.nzlmax = (std::max)(Index(4), fillratio) * annz / 4; // estimated nnz in L factor - // Return the estimated size to the user if necessary Index tempSpace; tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar); @@ -166,14 +178,10 @@ Index SparseLUImpl::memInit(Index m, Index n, Index annz, Index lw // Reserve memory for L/U factors do { - try - { - expand(glu.lusup, glu.nzlumax, 0, 0, num_expansions); - expand(glu.ucol,glu.nzumax, 0, 0, num_expansions); - expand(glu.lsub,glu.nzlmax, 0, 0, num_expansions); - expand(glu.usub,glu.nzumax, 0, 1, num_expansions); - } - catch(std::bad_alloc& ) + if( (expand(glu.lusup, glu.nzlumax, 0, 0, num_expansions)<0) + || (expand(glu.ucol, glu.nzumax, 0, 0, num_expansions)<0) + || (expand (glu.lsub, glu.nzlmax, 0, 0, num_expansions)<0) + || (expand (glu.usub, glu.nzumax, 0, 1, num_expansions)<0) ) { //Reduce the estimated size and retry glu.nzlumax /= 2; @@ -181,10 +189,7 @@ Index SparseLUImpl::memInit(Index m, Index n, Index annz, Index lw glu.nzlmax /= 2; if (glu.nzlumax < annz ) return glu.nzlumax; } - } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size()); - - ++num_expansions; return 0; diff --git a/otherlibs/eigen3/Eigen/src/SparseQR/SparseQR.h b/otherlibs/eigen3/Eigen/src/SparseQR/SparseQR.h index 07c46e4b90..afda43bfc6 100644 --- a/otherlibs/eigen3/Eigen/src/SparseQR/SparseQR.h +++ b/otherlibs/eigen3/Eigen/src/SparseQR/SparseQR.h @@ -161,8 +161,9 @@ class SparseQR b = y; // Solve with the triangular matrix R + y.resize((std::max)(cols(),Index(y.rows())),y.cols()); y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView().solve(b.topRows(rank)); - y.bottomRows(y.size()-rank).setZero(); + y.bottomRows(y.rows()-rank).setZero(); // Apply the column permutation if (m_perm_c.size()) dest.topRows(cols()) = colsPermutation() * y.topRows(cols()); @@ -246,7 +247,7 @@ class SparseQR Index m_nonzeropivots; // Number of non zero pivots found 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_isQSorted; // whether Q is sorted or not template friend struct SparseQR_QProduct; template friend struct SparseQRMatrixQReturnType; @@ -338,7 +339,7 @@ void SparseQR::factorize(const MatrixType& mat) Index nonzeroCol = 0; // Record the number of valid pivots // Left looking rank-revealing QR factorization: compute a column of R and Q at a time - for (Index col = 0; col < n; ++col) + for (Index col = 0; col < (std::min)(n,m); ++col) { mark.setConstant(-1); m_R.startVec(col); @@ -346,7 +347,7 @@ void SparseQR::factorize(const MatrixType& mat) mark(nonzeroCol) = col; Qidx(0) = nonzeroCol; nzcolR = 0; nzcolQ = 1; - found_diag = false; + found_diag = col>=m; tval.setZero(); // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., @@ -571,6 +572,7 @@ struct SparseQR_QProduct : ReturnByValue topRightCorner() const /** \returns an expression of a top-right corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -188,13 +188,13 @@ inline const Block topLeftCorner() const /** \returns an expression of a top-left corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -263,13 +263,13 @@ inline const Block bottomRightCorner() const /** \returns an expression of a bottom-right corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -338,13 +338,13 @@ inline const Block bottomLeftCorner() const /** \returns an expression of a bottom-left corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -390,7 +390,11 @@ inline ConstRowsBlockXpr topRows(Index n) const /** \returns a block consisting of the top rows of *this. * - * \tparam N the number of rows in the block + * \tparam N the number of rows in the block as specified at compile-time + * \param n the number of rows in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_topRows.cpp * Output: \verbinclude MatrixBase_template_int_topRows.out @@ -398,16 +402,16 @@ inline ConstRowsBlockXpr topRows(Index n) const * \sa class Block, block(Index,Index,Index,Index) */ template -inline typename NRowsBlockXpr::Type topRows() +inline typename NRowsBlockXpr::Type topRows(Index n = N) { - return typename NRowsBlockXpr::Type(derived(), 0, 0, N, cols()); + return typename NRowsBlockXpr::Type(derived(), 0, 0, n, cols()); } /** This is the const version of topRows().*/ template -inline typename ConstNRowsBlockXpr::Type topRows() const +inline typename ConstNRowsBlockXpr::Type topRows(Index n = N) const { - return typename ConstNRowsBlockXpr::Type(derived(), 0, 0, N, cols()); + return typename ConstNRowsBlockXpr::Type(derived(), 0, 0, n, cols()); } @@ -434,7 +438,11 @@ inline ConstRowsBlockXpr bottomRows(Index n) const /** \returns a block consisting of the bottom rows of *this. * - * \tparam N the number of rows in the block + * \tparam N the number of rows in the block as specified at compile-time + * \param n the number of rows in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_bottomRows.cpp * Output: \verbinclude MatrixBase_template_int_bottomRows.out @@ -442,16 +450,16 @@ inline ConstRowsBlockXpr bottomRows(Index n) const * \sa class Block, block(Index,Index,Index,Index) */ template -inline typename NRowsBlockXpr::Type bottomRows() +inline typename NRowsBlockXpr::Type bottomRows(Index n = N) { - return typename NRowsBlockXpr::Type(derived(), rows() - N, 0, N, cols()); + return typename NRowsBlockXpr::Type(derived(), rows() - n, 0, n, cols()); } /** This is the const version of bottomRows().*/ template -inline typename ConstNRowsBlockXpr::Type bottomRows() const +inline typename ConstNRowsBlockXpr::Type bottomRows(Index n = N) const { - return typename ConstNRowsBlockXpr::Type(derived(), rows() - N, 0, N, cols()); + return typename ConstNRowsBlockXpr::Type(derived(), rows() - n, 0, n, cols()); } @@ -459,28 +467,32 @@ inline typename ConstNRowsBlockXpr::Type bottomRows() const /** \returns a block consisting of a range of rows of *this. * * \param startRow the index of the first row in the block - * \param numRows the number of rows in the block + * \param n the number of rows in the block * * Example: \include DenseBase_middleRows_int.cpp * Output: \verbinclude DenseBase_middleRows_int.out * * \sa class Block, block(Index,Index,Index,Index) */ -inline RowsBlockXpr middleRows(Index startRow, Index numRows) +inline RowsBlockXpr middleRows(Index startRow, Index n) { - return RowsBlockXpr(derived(), startRow, 0, numRows, cols()); + return RowsBlockXpr(derived(), startRow, 0, n, cols()); } /** This is the const version of middleRows(Index,Index).*/ -inline ConstRowsBlockXpr middleRows(Index startRow, Index numRows) const +inline ConstRowsBlockXpr middleRows(Index startRow, Index n) const { - return ConstRowsBlockXpr(derived(), startRow, 0, numRows, cols()); + return ConstRowsBlockXpr(derived(), startRow, 0, n, cols()); } /** \returns a block consisting of a range of rows of *this. * - * \tparam N the number of rows in the block + * \tparam N the number of rows in the block as specified at compile-time * \param startRow the index of the first row in the block + * \param n the number of rows in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include DenseBase_template_int_middleRows.cpp * Output: \verbinclude DenseBase_template_int_middleRows.out @@ -488,16 +500,16 @@ inline ConstRowsBlockXpr middleRows(Index startRow, Index numRows) const * \sa class Block, block(Index,Index,Index,Index) */ template -inline typename NRowsBlockXpr::Type middleRows(Index startRow) +inline typename NRowsBlockXpr::Type middleRows(Index startRow, Index n = N) { - return typename NRowsBlockXpr::Type(derived(), startRow, 0, N, cols()); + return typename NRowsBlockXpr::Type(derived(), startRow, 0, n, cols()); } /** This is the const version of middleRows().*/ template -inline typename ConstNRowsBlockXpr::Type middleRows(Index startRow) const +inline typename ConstNRowsBlockXpr::Type middleRows(Index startRow, Index n = N) const { - return typename ConstNRowsBlockXpr::Type(derived(), startRow, 0, N, cols()); + return typename ConstNRowsBlockXpr::Type(derived(), startRow, 0, n, cols()); } @@ -524,7 +536,11 @@ inline ConstColsBlockXpr leftCols(Index n) const /** \returns a block consisting of the left columns of *this. * - * \tparam N the number of columns in the block + * \tparam N the number of columns in the block as specified at compile-time + * \param n the number of columns in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_leftCols.cpp * Output: \verbinclude MatrixBase_template_int_leftCols.out @@ -532,16 +548,16 @@ inline ConstColsBlockXpr leftCols(Index n) const * \sa class Block, block(Index,Index,Index,Index) */ template -inline typename NColsBlockXpr::Type leftCols() +inline typename NColsBlockXpr::Type leftCols(Index n = N) { - return typename NColsBlockXpr::Type(derived(), 0, 0, rows(), N); + return typename NColsBlockXpr::Type(derived(), 0, 0, rows(), n); } /** This is the const version of leftCols().*/ template -inline typename ConstNColsBlockXpr::Type leftCols() const +inline typename ConstNColsBlockXpr::Type leftCols(Index n = N) const { - return typename ConstNColsBlockXpr::Type(derived(), 0, 0, rows(), N); + return typename ConstNColsBlockXpr::Type(derived(), 0, 0, rows(), n); } @@ -568,7 +584,11 @@ inline ConstColsBlockXpr rightCols(Index n) const /** \returns a block consisting of the right columns of *this. * - * \tparam N the number of columns in the block + * \tparam N the number of columns in the block as specified at compile-time + * \param n the number of columns in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_rightCols.cpp * Output: \verbinclude MatrixBase_template_int_rightCols.out @@ -576,16 +596,16 @@ inline ConstColsBlockXpr rightCols(Index n) const * \sa class Block, block(Index,Index,Index,Index) */ template -inline typename NColsBlockXpr::Type rightCols() +inline typename NColsBlockXpr::Type rightCols(Index n = N) { - return typename NColsBlockXpr::Type(derived(), 0, cols() - N, rows(), N); + return typename NColsBlockXpr::Type(derived(), 0, cols() - n, rows(), n); } /** This is the const version of rightCols().*/ template -inline typename ConstNColsBlockXpr::Type rightCols() const +inline typename ConstNColsBlockXpr::Type rightCols(Index n = N) const { - return typename ConstNColsBlockXpr::Type(derived(), 0, cols() - N, rows(), N); + return typename ConstNColsBlockXpr::Type(derived(), 0, cols() - n, rows(), n); } @@ -613,8 +633,12 @@ inline ConstColsBlockXpr middleCols(Index startCol, Index numCols) const /** \returns a block consisting of a range of columns of *this. * - * \tparam N the number of columns in the block + * \tparam N the number of columns in the block as specified at compile-time * \param startCol the index of the first column in the block + * \param n the number of columns in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include DenseBase_template_int_middleCols.cpp * Output: \verbinclude DenseBase_template_int_middleCols.out @@ -622,16 +646,16 @@ inline ConstColsBlockXpr middleCols(Index startCol, Index numCols) const * \sa class Block, block(Index,Index,Index,Index) */ template -inline typename NColsBlockXpr::Type middleCols(Index startCol) +inline typename NColsBlockXpr::Type middleCols(Index startCol, Index n = N) { - return typename NColsBlockXpr::Type(derived(), 0, startCol, rows(), N); + return typename NColsBlockXpr::Type(derived(), 0, startCol, rows(), n); } /** This is the const version of middleCols().*/ template -inline typename ConstNColsBlockXpr::Type middleCols(Index startCol) const +inline typename ConstNColsBlockXpr::Type middleCols(Index startCol, Index n = N) const { - return typename ConstNColsBlockXpr::Type(derived(), 0, startCol, rows(), N); + return typename ConstNColsBlockXpr::Type(derived(), 0, startCol, rows(), n); } @@ -667,15 +691,15 @@ inline const Block block(Index startRow, In /** \returns an expression of a block in *this. * - * \tparam BlockRows number of rows in block as specified at compile time - * \tparam BlockCols number of columns in block as specified at compile time + * \tparam BlockRows number of rows in block as specified at compile-time + * \tparam BlockCols number of columns in block as specified at compile-time * \param startRow the first row in the block * \param startCol the first column in the block - * \param blockRows number of rows in block as specified at run time - * \param blockCols number of columns in block as specified at run time + * \param blockRows number of rows in block as specified at run-time + * \param blockCols number of columns in block as specified at run-time * - * This function is mainly useful for blocks where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for blocks where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a blockRows should equal \a BlockRows unless * \a BlockRows is \a Dynamic, and the same for the number of columns. * @@ -738,7 +762,7 @@ inline ConstRowXpr row(Index i) const * \only_for_vectors * * \param start the first coefficient in the segment - * \param vecSize the number of coefficients in the segment + * \param n the number of coefficients in the segment * * Example: \include MatrixBase_segment_int_int.cpp * Output: \verbinclude MatrixBase_segment_int_int.out @@ -749,25 +773,25 @@ inline ConstRowXpr row(Index i) const * * \sa class Block, segment(Index) */ -inline SegmentReturnType segment(Index start, Index vecSize) +inline SegmentReturnType segment(Index start, Index n) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return SegmentReturnType(derived(), start, vecSize); + return SegmentReturnType(derived(), start, n); } /** This is the const version of segment(Index,Index).*/ -inline ConstSegmentReturnType segment(Index start, Index vecSize) const +inline ConstSegmentReturnType segment(Index start, Index n) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return ConstSegmentReturnType(derived(), start, vecSize); + return ConstSegmentReturnType(derived(), start, n); } /** \returns a dynamic-size expression of the first coefficients of *this. * * \only_for_vectors * - * \param vecSize the number of coefficients in the block + * \param n the number of coefficients in the segment * * Example: \include MatrixBase_start_int.cpp * Output: \verbinclude MatrixBase_start_int.out @@ -778,25 +802,24 @@ inline ConstSegmentReturnType segment(Index start, Index vecSize) const * * \sa class Block, block(Index,Index) */ -inline SegmentReturnType head(Index vecSize) +inline SegmentReturnType head(Index n) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return SegmentReturnType(derived(), 0, vecSize); + return SegmentReturnType(derived(), 0, n); } /** This is the const version of head(Index).*/ -inline ConstSegmentReturnType - head(Index vecSize) const +inline ConstSegmentReturnType head(Index n) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return ConstSegmentReturnType(derived(), 0, vecSize); + return ConstSegmentReturnType(derived(), 0, n); } /** \returns a dynamic-size expression of the last coefficients of *this. * * \only_for_vectors * - * \param vecSize the number of coefficients in the block + * \param n the number of coefficients in the segment * * Example: \include MatrixBase_end_int.cpp * Output: \verbinclude MatrixBase_end_int.out @@ -807,95 +830,106 @@ inline ConstSegmentReturnType * * \sa class Block, block(Index,Index) */ -inline SegmentReturnType tail(Index vecSize) +inline SegmentReturnType tail(Index n) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return SegmentReturnType(derived(), this->size() - vecSize, vecSize); + return SegmentReturnType(derived(), this->size() - n, n); } /** This is the const version of tail(Index).*/ -inline ConstSegmentReturnType tail(Index vecSize) const +inline ConstSegmentReturnType tail(Index n) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return ConstSegmentReturnType(derived(), this->size() - vecSize, vecSize); + return ConstSegmentReturnType(derived(), this->size() - n, n); } /** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this * * \only_for_vectors * - * The template parameter \a Size is the number of coefficients in the block + * \tparam N the number of coefficients in the segment as specified at compile-time + * \param start the index of the first element in the segment + * \param n the number of coefficients in the segment as specified at compile-time * - * \param start the index of the first element of the sub-vector + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_segment.cpp * Output: \verbinclude MatrixBase_template_int_segment.out * * \sa class Block */ -template -inline typename FixedSegmentReturnType::Type segment(Index start) +template +inline typename FixedSegmentReturnType::Type segment(Index start, Index n = N) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename FixedSegmentReturnType::Type(derived(), start); + return typename FixedSegmentReturnType::Type(derived(), start, n); } /** This is the const version of segment(Index).*/ -template -inline typename ConstFixedSegmentReturnType::Type segment(Index start) const +template +inline typename ConstFixedSegmentReturnType::Type segment(Index start, Index n = N) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename ConstFixedSegmentReturnType::Type(derived(), start); + return typename ConstFixedSegmentReturnType::Type(derived(), start, n); } /** \returns a fixed-size expression of the first coefficients of *this. * * \only_for_vectors * - * The template parameter \a Size is the number of coefficients in the block + * \tparam N the number of coefficients in the segment as specified at compile-time + * \param n the number of coefficients in the segment as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_start.cpp * Output: \verbinclude MatrixBase_template_int_start.out * * \sa class Block */ -template -inline typename FixedSegmentReturnType::Type head() +template +inline typename FixedSegmentReturnType::Type head(Index n = N) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename FixedSegmentReturnType::Type(derived(), 0); + return typename FixedSegmentReturnType::Type(derived(), 0, n); } /** This is the const version of head().*/ -template -inline typename ConstFixedSegmentReturnType::Type head() const +template +inline typename ConstFixedSegmentReturnType::Type head(Index n = N) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename ConstFixedSegmentReturnType::Type(derived(), 0); + return typename ConstFixedSegmentReturnType::Type(derived(), 0, n); } /** \returns a fixed-size expression of the last coefficients of *this. * * \only_for_vectors * - * The template parameter \a Size is the number of coefficients in the block + * \tparam N the number of coefficients in the segment as specified at compile-time + * \param n the number of coefficients in the segment as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_end.cpp * Output: \verbinclude MatrixBase_template_int_end.out * * \sa class Block */ -template -inline typename FixedSegmentReturnType::Type tail() +template +inline typename FixedSegmentReturnType::Type tail(Index n = N) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename FixedSegmentReturnType::Type(derived(), size() - Size); + return typename FixedSegmentReturnType::Type(derived(), size() - n); } /** This is the const version of tail.*/ -template -inline typename ConstFixedSegmentReturnType::Type tail() const +template +inline typename ConstFixedSegmentReturnType::Type tail(Index n = N) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename ConstFixedSegmentReturnType::Type(derived(), size() - Size); + return typename ConstFixedSegmentReturnType::Type(derived(), size() - n); } diff --git a/otherlibs/eigen3/Eigen/src/plugins/MatrixCwiseBinaryOps.h b/otherlibs/eigen3/Eigen/src/plugins/MatrixCwiseBinaryOps.h index 3a737df7b8..7f62149e04 100644 --- a/otherlibs/eigen3/Eigen/src/plugins/MatrixCwiseBinaryOps.h +++ b/otherlibs/eigen3/Eigen/src/plugins/MatrixCwiseBinaryOps.h @@ -83,7 +83,7 @@ cwiseMin(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const EIGEN_STRONG_INLINE const CwiseBinaryOp, const Derived, const ConstantReturnType> cwiseMin(const Scalar &other) const { - return cwiseMin(Derived::PlainObject::Constant(rows(), cols(), other)); + return cwiseMin(Derived::Constant(rows(), cols(), other)); } /** \returns an expression of the coefficient-wise max of *this and \a other @@ -107,7 +107,7 @@ cwiseMax(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const EIGEN_STRONG_INLINE const CwiseBinaryOp, const Derived, const ConstantReturnType> cwiseMax(const Scalar &other) const { - return cwiseMax(Derived::PlainObject::Constant(rows(), cols(), other)); + return cwiseMax(Derived::Constant(rows(), cols(), other)); } diff --git a/otherlibs/eigen3/unsupported/Eigen/AdolcForward b/otherlibs/eigen3/unsupported/Eigen/AdolcForward index e40b3af47e..2627decd0f 100644 --- a/otherlibs/eigen3/unsupported/Eigen/AdolcForward +++ b/otherlibs/eigen3/unsupported/Eigen/AdolcForward @@ -45,8 +45,7 @@ namespace Eigen { /** - * \defgroup AdolcForward_Module Adolc forward module - * \ingroup eigen_grp + * \defgroup AdolcForward_Module Adolc forward module * This module provides support for adolc's adouble type in forward mode. * ADOL-C is a C++ automatic differentiation library, * see https://projects.coin-or.org/ADOL-C for more information. diff --git a/otherlibs/eigen3/unsupported/Eigen/AlignedVector3 b/otherlibs/eigen3/unsupported/Eigen/AlignedVector3 index 9d33aabd3b..7b45e6cce1 100644 --- a/otherlibs/eigen3/unsupported/Eigen/AlignedVector3 +++ b/otherlibs/eigen3/unsupported/Eigen/AlignedVector3 @@ -15,8 +15,7 @@ namespace Eigen { /** - * \defgroup AlignedVector3_Module Aligned vector3 module - * \ingroup eigen_grp + * \defgroup AlignedVector3_Module Aligned vector3 module * * \code * #include diff --git a/otherlibs/eigen3/unsupported/Eigen/ArpackSupport b/otherlibs/eigen3/unsupported/Eigen/ArpackSupport index f0cc70cb1b..37a2799ef2 100644 --- a/otherlibs/eigen3/unsupported/Eigen/ArpackSupport +++ b/otherlibs/eigen3/unsupported/Eigen/ArpackSupport @@ -13,8 +13,7 @@ #include -/** \defgroup ArpackSupport_Module Arpack support module - * \ingroup eigen_grp +/** \defgroup ArpackSupport_Module Arpack support module * * This module provides a wrapper to Arpack, a library for sparse eigenvalue decomposition. * diff --git a/otherlibs/eigen3/unsupported/Eigen/AutoDiff b/otherlibs/eigen3/unsupported/Eigen/AutoDiff index 27ed8c497a..abf5b7d678 100644 --- a/otherlibs/eigen3/unsupported/Eigen/AutoDiff +++ b/otherlibs/eigen3/unsupported/Eigen/AutoDiff @@ -13,8 +13,7 @@ namespace Eigen { /** - * \defgroup AutoDiff_Module Auto Diff module - * \ingroup eigen_grp + * \defgroup AutoDiff_Module Auto Diff module * * This module features forward automatic differentation via a simple * templated scalar type wrapper AutoDiffScalar. diff --git a/otherlibs/eigen3/unsupported/Eigen/BVH b/otherlibs/eigen3/unsupported/Eigen/BVH index 864ed99349..0161a5402d 100644 --- a/otherlibs/eigen3/unsupported/Eigen/BVH +++ b/otherlibs/eigen3/unsupported/Eigen/BVH @@ -19,8 +19,7 @@ namespace Eigen { /** - * \defgroup BVH_Module BVH module - * \ingroup eigen_grp + * \defgroup BVH_Module BVH module * \brief This module provides generic bounding volume hierarchy algorithms * and reference tree implementations. * diff --git a/otherlibs/eigen3/unsupported/Eigen/FFT b/otherlibs/eigen3/unsupported/Eigen/FFT index 3b7fa23b62..2c45b3999e 100644 --- a/otherlibs/eigen3/unsupported/Eigen/FFT +++ b/otherlibs/eigen3/unsupported/Eigen/FFT @@ -17,8 +17,7 @@ /** - * \defgroup FFT_Module Fast Fourier Transform module - * \ingroup eigen_grp + * \defgroup FFT_Module Fast Fourier Transform module * * \code * #include diff --git a/otherlibs/eigen3/unsupported/Eigen/IterativeSolvers b/otherlibs/eigen3/unsupported/Eigen/IterativeSolvers index 91b37f01ef..aa15403dbf 100644 --- a/otherlibs/eigen3/unsupported/Eigen/IterativeSolvers +++ b/otherlibs/eigen3/unsupported/Eigen/IterativeSolvers @@ -13,8 +13,7 @@ #include /** - * \defgroup IterativeSolvers_Module Iterative solvers module - * \ingroup eigen_grp + * \defgroup IterativeSolvers_Module Iterative solvers module * This module aims to provide various iterative linear and non linear solver algorithms. * It currently provides: * - a constrained conjugate gradient diff --git a/otherlibs/eigen3/unsupported/Eigen/KroneckerProduct b/otherlibs/eigen3/unsupported/Eigen/KroneckerProduct index 10e5840155..c932c06a6d 100644 --- a/otherlibs/eigen3/unsupported/Eigen/KroneckerProduct +++ b/otherlibs/eigen3/unsupported/Eigen/KroneckerProduct @@ -16,8 +16,7 @@ namespace Eigen { /** - * \defgroup KroneckerProduct_Module KroneckerProduct module - * \ingroup eigen_grp + * \defgroup KroneckerProduct_Module KroneckerProduct module * * This module contains an experimental Kronecker product implementation. * diff --git a/otherlibs/eigen3/unsupported/Eigen/LevenbergMarquardt b/otherlibs/eigen3/unsupported/Eigen/LevenbergMarquardt index af4dc5fe13..0fe2680bab 100644 --- a/otherlibs/eigen3/unsupported/Eigen/LevenbergMarquardt +++ b/otherlibs/eigen3/unsupported/Eigen/LevenbergMarquardt @@ -20,8 +20,7 @@ #include /** - * \defgroup LevenbergMarquardt_Module Levenberg-Marquardt module - * \ingroup eigen_grp + * \defgroup LevenbergMarquardt_Module Levenberg-Marquardt module * * \code * #include diff --git a/otherlibs/eigen3/unsupported/Eigen/MPRealSupport b/otherlibs/eigen3/unsupported/Eigen/MPRealSupport index 4eb96ad81e..d4b03647db 100644 --- a/otherlibs/eigen3/unsupported/Eigen/MPRealSupport +++ b/otherlibs/eigen3/unsupported/Eigen/MPRealSupport @@ -18,8 +18,7 @@ namespace Eigen { /** - * \defgroup MPRealSupport_Module MPFRC++ Support module - * \ingroup eigen_grp + * \defgroup MPRealSupport_Module MPFRC++ Support module * \code * #include * \endcode diff --git a/otherlibs/eigen3/unsupported/Eigen/MatrixFunctions b/otherlibs/eigen3/unsupported/Eigen/MatrixFunctions index 8ef2d1d47b..0991817d55 100644 --- a/otherlibs/eigen3/unsupported/Eigen/MatrixFunctions +++ b/otherlibs/eigen3/unsupported/Eigen/MatrixFunctions @@ -21,8 +21,7 @@ #include /** - * \defgroup MatrixFunctions_Module Matrix functions module - * \ingroup eigen_grp + * \defgroup MatrixFunctions_Module Matrix functions module * \brief This module aims to provide various methods for the computation of * matrix functions. * diff --git a/otherlibs/eigen3/unsupported/Eigen/MoreVectorization b/otherlibs/eigen3/unsupported/Eigen/MoreVectorization index dbc5750f2b..470e724308 100644 --- a/otherlibs/eigen3/unsupported/Eigen/MoreVectorization +++ b/otherlibs/eigen3/unsupported/Eigen/MoreVectorization @@ -14,8 +14,7 @@ namespace Eigen { /** - * \defgroup MoreVectorization More vectorization module - * \ingroup eigen_grp + * \defgroup MoreVectorization More vectorization module */ } diff --git a/otherlibs/eigen3/unsupported/Eigen/NonLinearOptimization b/otherlibs/eigen3/unsupported/Eigen/NonLinearOptimization index 0005cba058..600ab4c12e 100644 --- a/otherlibs/eigen3/unsupported/Eigen/NonLinearOptimization +++ b/otherlibs/eigen3/unsupported/Eigen/NonLinearOptimization @@ -18,8 +18,7 @@ #include /** - * \defgroup NonLinearOptimization_Module Non linear optimization module - * \ingroup eigen_grp + * \defgroup NonLinearOptimization_Module Non linear optimization module * * \code * #include diff --git a/otherlibs/eigen3/unsupported/Eigen/NumericalDiff b/otherlibs/eigen3/unsupported/Eigen/NumericalDiff index e1776d31e3..433334ca80 100644 --- a/otherlibs/eigen3/unsupported/Eigen/NumericalDiff +++ b/otherlibs/eigen3/unsupported/Eigen/NumericalDiff @@ -15,8 +15,7 @@ namespace Eigen { /** - * \defgroup NumericalDiff_Module Numerical differentiation module - * \ingroup eigen_grp + * \defgroup NumericalDiff_Module Numerical differentiation module * * \code * #include diff --git a/otherlibs/eigen3/unsupported/Eigen/OpenGLSupport b/otherlibs/eigen3/unsupported/Eigen/OpenGLSupport index a6d5565c1d..c4090ab11d 100644 --- a/otherlibs/eigen3/unsupported/Eigen/OpenGLSupport +++ b/otherlibs/eigen3/unsupported/Eigen/OpenGLSupport @@ -11,13 +11,17 @@ #define EIGEN_OPENGL_MODULE #include -#include + +#if defined(__APPLE_CC__) + #include +#else + #include +#endif namespace Eigen { /** - * \defgroup OpenGLSUpport_Module OpenGL Support module - * \ingroup eigen_grp + * \defgroup OpenGLSUpport_Module OpenGL Support module * * This module provides wrapper functions for a couple of OpenGL functions * which simplify the way to pass Eigen's object to openGL. diff --git a/otherlibs/eigen3/unsupported/Eigen/Polynomials b/otherlibs/eigen3/unsupported/Eigen/Polynomials index ddda12e0a1..cece563374 100644 --- a/otherlibs/eigen3/unsupported/Eigen/Polynomials +++ b/otherlibs/eigen3/unsupported/Eigen/Polynomials @@ -25,8 +25,7 @@ #endif /** - * \defgroup Polynomials_Module Polynomials module - * \ingroup eigen_grp + * \defgroup Polynomials_Module Polynomials module * \brief This module provides a QR based polynomial solver. * * To use this module, add diff --git a/otherlibs/eigen3/unsupported/Eigen/SVD b/otherlibs/eigen3/unsupported/Eigen/SVD index 2e1dd608d7..7cc0592805 100644 --- a/otherlibs/eigen3/unsupported/Eigen/SVD +++ b/otherlibs/eigen3/unsupported/Eigen/SVD @@ -7,8 +7,7 @@ #include "../../Eigen/src/Core/util/DisableStupidWarnings.h" -/** \defgroup SVD_Module SVD module - * \ingroup eigen_grp +/** \defgroup SVD_Module SVD module * * * diff --git a/otherlibs/eigen3/unsupported/Eigen/Skyline b/otherlibs/eigen3/unsupported/Eigen/Skyline index aa255a2eb1..71a68cb422 100644 --- a/otherlibs/eigen3/unsupported/Eigen/Skyline +++ b/otherlibs/eigen3/unsupported/Eigen/Skyline @@ -20,8 +20,7 @@ #include /** - * \defgroup Skyline_Module Skyline module - * \ingroup eigen_grp + * \defgroup Skyline_Module Skyline module * * * diff --git a/otherlibs/eigen3/unsupported/Eigen/SparseExtra b/otherlibs/eigen3/unsupported/Eigen/SparseExtra index ccc6aaca8a..b5597902af 100644 --- a/otherlibs/eigen3/unsupported/Eigen/SparseExtra +++ b/otherlibs/eigen3/unsupported/Eigen/SparseExtra @@ -27,8 +27,7 @@ #endif /** - * \defgroup SparseExtra_Module SparseExtra module - * \ingroup eigen_grp + * \defgroup SparseExtra_Module SparseExtra module * * This module contains some experimental features extending the sparse module. * diff --git a/otherlibs/eigen3/unsupported/Eigen/Splines b/otherlibs/eigen3/unsupported/Eigen/Splines index 59c7716045..322e6b9f5c 100644 --- a/otherlibs/eigen3/unsupported/Eigen/Splines +++ b/otherlibs/eigen3/unsupported/Eigen/Splines @@ -13,8 +13,7 @@ namespace Eigen { /** - * \defgroup Splines_Module Spline and spline fitting module - * \ingroup eigen_grp + * \defgroup Splines_Module Spline and spline fitting module * * This module provides a simple multi-dimensional spline class while * offering most basic functionality to fit a spline to point sets. diff --git a/otherlibs/eigen3/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h b/otherlibs/eigen3/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h index 3f18beeebe..dc0093eb97 100644 --- a/otherlibs/eigen3/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h +++ b/otherlibs/eigen3/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h @@ -58,7 +58,9 @@ void pseudo_inverse(const CMatrix &C, CINVMatrix &CINV) Scalar rho, rho_1, alpha; d.setZero(); - CINV.startFill(); // FIXME estimate the number of non-zeros + typedef Triplet T; + std::vector tripletList; + for (Index i = 0; i < rows; ++i) { d[i] = 1.0; @@ -84,11 +86,12 @@ void pseudo_inverse(const CMatrix &C, CINVMatrix &CINV) // FIXME add a generic "prune/filter" expression for both dense and sparse object to sparse for (Index j=0; j TmpVec; diff --git a/otherlibs/eigen3/unsupported/doc/Overview.dox b/otherlibs/eigen3/unsupported/doc/Overview.dox index 0b3a95ee06..d048377df4 100644 --- a/otherlibs/eigen3/unsupported/doc/Overview.dox +++ b/otherlibs/eigen3/unsupported/doc/Overview.dox @@ -14,8 +14,7 @@ Don't miss the official Eigen documentation. /* -\defgroup Unsupported_modules Unsupported modules - * \ingroup eigen_grp +\defgroup Unsupported_modules Unsupported modules The unsupported modules are contributions from various users. They are provided "as is", without any support. Nevertheless, some of them are diff --git a/otherlibs/eigen3/unsupported/test/FFTW.cpp b/otherlibs/eigen3/unsupported/test/FFTW.cpp index a07bf274b8..d3718e2d22 100644 --- a/otherlibs/eigen3/unsupported/test/FFTW.cpp +++ b/otherlibs/eigen3/unsupported/test/FFTW.cpp @@ -16,9 +16,6 @@ std::complex RandomCpx() { return std::complex( (T)(rand()/(T)RAND_MAX - . using namespace std; using namespace Eigen; -float norm(float x) {return x*x;} -double norm(double x) {return x*x;} -long double norm(long double x) {return x*x;} template < typename T> complex promote(complex x) { return complex(x.real(),x.imag()); } @@ -40,11 +37,11 @@ complex promote(long double x) { return complex( x); for (size_t k1=0;k1<(size_t)timebuf.size();++k1) { acc += promote( timebuf[k1] ) * exp( complex(0,k1*phinc) ); } - totalpower += norm(acc); + totalpower += numext::abs2(acc); complex x = promote(fftbuf[k0]); complex dif = acc - x; - difpower += norm(dif); - //cerr << k0 << "\t" << acc << "\t" << x << "\t" << sqrt(norm(dif)) << endl; + difpower += numext::abs2(dif); + //cerr << k0 << "\t" << acc << "\t" << x << "\t" << sqrt(numext::abs2(dif)) << endl; } cerr << "rmse:" << sqrt(difpower/totalpower) << endl; return sqrt(difpower/totalpower); @@ -57,8 +54,8 @@ complex promote(long double x) { return complex( x); long double difpower=0; size_t n = (min)( buf1.size(),buf2.size() ); for (size_t k=0;k