From 5f8d58f36a57b860d6af52574df6f3e62debe001 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 9 Aug 2009 00:05:38 +0200 Subject: [PATCH 01/58] fix bug in sorting of singular values --- Eigen/src/SVD/SVD.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 9b7d955c7..acc567f94 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -393,8 +393,9 @@ void SVD::compute(const MatrixType& matrix) { int k; W.end(n-i).minCoeff(&k); - if (k != i) + if (k != 0) { + k += i; std::swap(W[k],W[i]); A.col(i).swap(A.col(k)); V.col(i).swap(V.col(k)); From fe813911f2fc3b5679c769a3624349a3e053ac63 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 9 Aug 2009 00:06:53 +0200 Subject: [PATCH 02/58] make LU::solve() not to crash when rank=0 --- Eigen/src/LU/LU.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index 43154b769..f36951fda 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -561,7 +561,7 @@ bool LU::solve( if(!isSurjective()) { // is c is in the image of U ? - RealScalar biggest_in_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : 0; for(int col = 0; col < c.cols(); ++col) for(int row = m_rank; row < c.rows(); ++row) if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision)) From 3ed83fa6813bd1401adce03a87ad452661d72f5e Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 9 Aug 2009 16:58:13 +0200 Subject: [PATCH 03/58] * add Jacobi transformations * add Jacobi (Hestenes) SVD decomposition for square matrices * add function for trivial Householder --- Eigen/Jacobi | 24 ++++ Eigen/SVD | 4 + Eigen/src/Core/MatrixBase.h | 6 + Eigen/src/Householder/Householder.h | 9 ++ Eigen/src/Jacobi/Jacobi.h | 91 +++++++++++++++ Eigen/src/SVD/JacobiSquareSVD.h | 169 ++++++++++++++++++++++++++++ 6 files changed, 303 insertions(+) create mode 100644 Eigen/Jacobi create mode 100644 Eigen/src/Jacobi/Jacobi.h create mode 100644 Eigen/src/SVD/JacobiSquareSVD.h diff --git a/Eigen/Jacobi b/Eigen/Jacobi new file mode 100644 index 000000000..33a6b757e --- /dev/null +++ b/Eigen/Jacobi @@ -0,0 +1,24 @@ +#ifndef EIGEN_JACOBI_MODULE_H +#define EIGEN_JACOBI_MODULE_H + +#include "Core" + +#include "src/Core/util/DisableMSVCWarnings.h" + +namespace Eigen { + +/** \defgroup Jacobi_Module Jacobi module + * This module provides Jacobi rotations. + * + * \code + * #include + * \endcode + */ + +#include "src/Jacobi/Jacobi.h" + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_JACOBI_MODULE_H diff --git a/Eigen/SVD b/Eigen/SVD index eef05564b..3aab35118 100644 --- a/Eigen/SVD +++ b/Eigen/SVD @@ -2,6 +2,8 @@ #define EIGEN_SVD_MODULE_H #include "Core" +#include "Householder" +#include "Jacobi" #include "src/Core/util/DisableMSVCWarnings.h" @@ -20,7 +22,9 @@ namespace Eigen { * \endcode */ +#include "src/SVD/Bidiagonalization.h" #include "src/SVD/SVD.h" +#include "src/SVD/JacobiSquareSVD.h" } // namespace Eigen diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 6ec7ddbb7..688b7c4d0 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -778,6 +778,12 @@ template class MatrixBase void applyHouseholderOnTheRight(const EssentialPart& essential, const RealScalar& beta); +///////// Jacobi module ///////// + + void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s); + void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s); + bool makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s); + bool makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s); #ifdef EIGEN_MATRIXBASE_PLUGIN #include EIGEN_MATRIXBASE_PLUGIN diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 8972806de..2a4b1395c 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -32,6 +32,15 @@ template struct ei_decrement_size }; }; +template +void makeTrivialHouseholder( + EssentialPart *essential, + typename EssentialPart::RealScalar *beta) +{ + *beta = typename EssentialPart::RealScalar(0); + essential->setZero(); +} + template template void MatrixBase::makeHouseholder( diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h new file mode 100644 index 000000000..993a723ab --- /dev/null +++ b/Eigen/src/Jacobi/Jacobi.h @@ -0,0 +1,91 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_JACOBI_H +#define EIGEN_JACOBI_H + +template +void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) +{ + for(int i = 0; i < cols(); ++i) + { + Scalar tmp = coeff(p,i); + coeffRef(p,i) = c * tmp - s * coeff(q,i); + coeffRef(q,i) = s * tmp + c * coeff(q,i); + } +} + +template +void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) +{ + for(int i = 0; i < rows(); ++i) + { + Scalar tmp = coeff(i,p); + coeffRef(i,p) = c * tmp - s * coeff(i,q); + coeffRef(i,q) = s * tmp + c * coeff(i,q); + } +} + +template +bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar max_coeff, Scalar *c, Scalar *s) +{ + if(ei_abs(y) < max_coeff * 0.5 * machine_epsilon()) + { + *c = Scalar(1); + *s = Scalar(0); + return true; + } + else + { + Scalar tau = (z - x) / (2 * y); + Scalar w = ei_sqrt(1 + ei_abs2(tau)); + Scalar t; + if(tau>0) + t = Scalar(1) / (tau + w); + else + t = Scalar(1) / (tau - w); + *c = Scalar(1) / ei_sqrt(1 + ei_abs2(t)); + *s = *c * t; + return false; + } +} + +template +inline bool MatrixBase::makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s) +{ + return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), max_coeff, c, s); +} + +template +inline bool MatrixBase::makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s) +{ + return ei_makeJacobi(col(p).squaredNorm(), + col(p).dot(col(q)), + col(q).squaredNorm(), + max_coeff, + c,s); +} + + +#endif // EIGEN_JACOBI_H diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h new file mode 100644 index 000000000..e9ecc89ac --- /dev/null +++ b/Eigen/src/SVD/JacobiSquareSVD.h @@ -0,0 +1,169 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Benoit Jacob +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_JACOBISQUARESVD_H +#define EIGEN_JACOBISQUARESVD_H + +/** \ingroup SVD_Module + * \nonstableyet + * + * \class JacobiSquareSVD + * + * \brief Jacobi SVD decomposition of a square matrix + * + * \param MatrixType the type of the matrix of which we are computing the SVD decomposition + * \param ComputeU whether the U matrix should be computed + * \param ComputeV whether the V matrix should be computed + * + * \sa MatrixBase::jacobiSvd() + */ +template class JacobiSquareSVD +{ + private: + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + Options = MatrixType::Options + }; + + typedef Matrix DummyMatrixType; + typedef typename ei_meta_if, + DummyMatrixType>::ret MatrixUType; + typedef typename Diagonal::PlainMatrixType SingularValuesType; + typedef Matrix RowType; + typedef Matrix ColType; + + public: + + JacobiSquareSVD() : m_isInitialized(false) {} + + JacobiSquareSVD(const MatrixType& matrix) + { + compute(matrix); + } + + void compute(const MatrixType& matrix); + + const MatrixUType& matrixU() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_matrixU; + } + + const SingularValuesType& singularValues() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_singularValues; + } + + const MatrixUType& matrixV() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_matrixV; + } + + protected: + MatrixUType m_matrixU; + MatrixUType m_matrixV; + SingularValuesType m_singularValues; + bool m_isInitialized; +}; + +template +void JacobiSquareSVD::compute(const MatrixType& matrix) +{ + MatrixType work_matrix(matrix); + int size = matrix.rows(); + if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); + if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); + m_singularValues.resize(size); + RealScalar max_coeff = work_matrix.cwise().abs().maxCoeff(); + for(int k = 1; k < 40; ++k) { + bool finished = true; + for(int p = 1; p < size; ++p) + { + for(int q = 0; q < p; ++q) + { + Scalar c, s; + finished &= work_matrix.makeJacobiForAtA(p,q,max_coeff,&c,&s); + work_matrix.applyJacobiOnTheRight(p,q,c,s); + if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s); + } + } + if(finished) break; + } + + for(int i = 0; i < size; ++i) + { + m_singularValues.coeffRef(i) = work_matrix.col(i).norm(); + } + + int first_zero = size; + RealScalar biggest = m_singularValues.maxCoeff(); + for(int i = 0; i < size; i++) + { + int pos; + RealScalar biggest_remaining = m_singularValues.end(size-i).maxCoeff(&pos); + if(first_zero == size && ei_isMuchSmallerThan(biggest_remaining, biggest)) first_zero = pos + i; + if(pos) + { + pos += i; + std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); + if(ComputeU) work_matrix.col(pos).swap(work_matrix.col(i)); + if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); + } + } + + if(ComputeU) + { + for(int i = 0; i < first_zero; ++i) + { + m_matrixU.col(i) = work_matrix.col(i) / m_singularValues.coeff(i); + } + if(first_zero < size) + { + for(int i = first_zero; i < size; ++i) + { + for(int j = 0; j < size; ++j) + { + m_matrixU.col(i).setZero(); + m_matrixU.coeffRef(j,i) = Scalar(1); + for(int k = 0; k < first_zero; ++k) + m_matrixU.col(i) -= m_matrixU.col(i).dot(m_matrixU.col(k)) * m_matrixU.col(k); + RealScalar n = m_matrixU.col(i).norm(); + if(!ei_isMuchSmallerThan(n, biggest)) + { + m_matrixU.col(i) /= n; + break; + } + } + } + } + } +} +#endif // EIGEN_JACOBISQUARESVD_H From 8e086801196a251c013479de51e6cdaa5f595a9b Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 9 Aug 2009 17:18:42 +0200 Subject: [PATCH 04/58] don't depend on uninitialized value --- Eigen/src/SVD/JacobiSquareSVD.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h index e9ecc89ac..ad55735fc 100644 --- a/Eigen/src/SVD/JacobiSquareSVD.h +++ b/Eigen/src/SVD/JacobiSquareSVD.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2009 Benoit Jacob +// Copyright (C) 2009 Benoit Jacob // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -62,7 +62,7 @@ template class JacobiSquareSV JacobiSquareSVD() : m_isInitialized(false) {} - JacobiSquareSVD(const MatrixType& matrix) + JacobiSquareSVD(const MatrixType& matrix) : m_isInitialized(false) { compute(matrix); } @@ -165,5 +165,6 @@ void JacobiSquareSVD::compute(const MatrixType& } } } + m_isInitialized = true; } #endif // EIGEN_JACOBISQUARESVD_H From 527557672ab11a99ceecc2a3bfd9fea9dbe216dd Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 9 Aug 2009 20:49:55 +0200 Subject: [PATCH 05/58] disable the assembly for fast unaligned stores. indeed, there is a strange bug that is triggered by this code: #include int main() { Eigen::Matrix4f m; m <<1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16; m.col(0).swap(m.col(1)); std::cout << m << std::endl; } when the fast unaligned stores are used, the column is copied instead of being swapped. --- Eigen/src/Core/arch/SSE/PacketMath.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 3fd33afbf..130aa63ce 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -178,7 +178,7 @@ template<> EIGEN_STRONG_INLINE Packet4f ei_pload(const float* from) { template<> EIGEN_STRONG_INLINE Packet2d ei_pload(const double* from) { return _mm_load_pd(from); } template<> EIGEN_STRONG_INLINE Packet4i ei_pload(const int* from) { return _mm_load_si128(reinterpret_cast(from)); } -#if (!defined __GNUC__) && (!defined __ICC) +#if 1 // (!defined __GNUC__) && (!defined __ICC) template<> EIGEN_STRONG_INLINE Packet4f ei_ploadu(const float* from) { return _mm_loadu_ps(from); } template<> EIGEN_STRONG_INLINE Packet2d ei_ploadu(const double* from) { return _mm_loadu_pd(from); } template<> EIGEN_STRONG_INLINE Packet4i ei_ploadu(const int* from) { return _mm_loadu_si128(reinterpret_cast(from)); } From 1f1705868b85ff9153beabb444c1061d3ad0b2d0 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 9 Aug 2009 21:35:13 +0200 Subject: [PATCH 06/58] now you can #define EIGEN_DEBUG_ASSIGN, and all the values in ei_assign_traits are printed --- Eigen/src/Core/Assign.h | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index f3521d3dd..c47280dff 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -90,6 +90,25 @@ public: ? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) ) : int(NoUnrolling) }; + + static void debug() + { + EIGEN_DEBUG_VAR(DstIsAligned) + EIGEN_DEBUG_VAR(SrcIsAligned) + EIGEN_DEBUG_VAR(SrcAlignment) + EIGEN_DEBUG_VAR(InnerSize) + EIGEN_DEBUG_VAR(InnerMaxSize) + EIGEN_DEBUG_VAR(PacketSize) + EIGEN_DEBUG_VAR(MightVectorize) + EIGEN_DEBUG_VAR(MayInnerVectorize) + EIGEN_DEBUG_VAR(MayLinearVectorize) + EIGEN_DEBUG_VAR(MaySliceVectorize) + EIGEN_DEBUG_VAR(Vectorization) + EIGEN_DEBUG_VAR(UnrollingLimit) + EIGEN_DEBUG_VAR(MayUnrollCompletely) + EIGEN_DEBUG_VAR(MayUnrollInner) + EIGEN_DEBUG_VAR(Unrolling) + } }; /*************************************************************************** @@ -405,6 +424,9 @@ EIGEN_STRONG_INLINE Derived& MatrixBase YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) ei_assert(rows() == other.rows() && cols() == other.cols()); ei_assign_impl::run(derived(),other.derived()); +#ifdef EIGEN_DEBUG_ASSIGN + ei_assign_traits::debug(); +#endif return derived(); } From 216ee335acb23cf27604c8a2a76013e897b96a05 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 9 Aug 2009 22:19:12 +0200 Subject: [PATCH 07/58] LinearVectorization: If the destination isn't aligned, we have to do runtime checks and we don't unroll, so it's only good for large enough sizes --- Eigen/src/Core/Assign.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index c47280dff..68b6c817d 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -57,7 +57,10 @@ private: && ((int(Derived::Flags)&RowMajorBit)==(int(OtherDerived::Flags)&RowMajorBit)), MayInnerVectorize = MightVectorize && int(InnerSize)!=Dynamic && int(InnerSize)%int(PacketSize)==0 && int(DstIsAligned) && int(SrcIsAligned), - MayLinearVectorize = MightVectorize && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit), + MayLinearVectorize = MightVectorize && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit) + && (DstIsAligned || InnerMaxSize == Dynamic),/* If the destination isn't aligned, + we have to do runtime checks and we don't unroll, so it's only good for large enough sizes. See remark below + about InnerMaxSize. */ MaySliceVectorize = MightVectorize && int(InnerMaxSize)>=3*PacketSize /* slice vectorization can be slow, so we only want it if the slices are big, which is indicated by InnerMaxSize rather than InnerSize, think of the case of a dynamic block in a fixed-size matrix */ From ef55e7f4cecba2b5157681eb1d518a59deb4c599 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 9 Aug 2009 23:09:46 +0200 Subject: [PATCH 08/58] make custom asm directive volatile --- Eigen/src/Core/arch/SSE/PacketMath.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 3fd33afbf..52e666d2f 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -191,22 +191,22 @@ template<> EIGEN_STRONG_INLINE Packet4i ei_ploadu(const int* from) { return template<> EIGEN_STRONG_INLINE Packet4f ei_ploadu(const float* from) { __m128 res; - asm("movsd %[from0], %[r]" : [r] "=x" (res) : [from0] "m" (*from), [dummy] "m" (*(from+1)) ); - asm("movhps %[from2], %[r]" : [r] "+x" (res) : [from2] "m" (*(from+2)), [dummy] "m" (*(from+3)) ); + asm volatile ("movsd %[from0], %[r]" : [r] "=x" (res) : [from0] "m" (*from), [dummy] "m" (*(from+1)) ); + asm volatile ("movhps %[from2], %[r]" : [r] "+x" (res) : [from2] "m" (*(from+2)), [dummy] "m" (*(from+3)) ); return res; } template<> EIGEN_STRONG_INLINE Packet2d ei_ploadu(const double* from) { __m128d res; - asm("movsd %[from0], %[r]" : [r] "=x" (res) : [from0] "m" (*from) ); - asm("movhpd %[from1], %[r]" : [r] "+x" (res) : [from1] "m" (*(from+1)) ); + asm volatile ("movsd %[from0], %[r]" : [r] "=x" (res) : [from0] "m" (*from) ); + asm volatile ("movhpd %[from1], %[r]" : [r] "+x" (res) : [from1] "m" (*(from+1)) ); return res; } template<> EIGEN_STRONG_INLINE Packet4i ei_ploadu(const int* from) { __m128i res; - asm("movsd %[from0], %[r]" : [r] "=x" (res) : [from0] "m" (*from), [dummy] "m" (*(from+1)) ); - asm("movhps %[from2], %[r]" : [r] "+x" (res) : [from2] "m" (*(from+2)), [dummy] "m" (*(from+3)) ); + asm volatile ("movsd %[from0], %[r]" : [r] "=x" (res) : [from0] "m" (*from), [dummy] "m" (*(from+1)) ); + asm volatile ("movhps %[from2], %[r]" : [r] "+x" (res) : [from2] "m" (*(from+2)), [dummy] "m" (*(from+3)) ); return res; } #endif From ea884e6f48380c4ec8e26f8cb9af98c8f9d12db6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 11 Aug 2009 15:08:03 +0200 Subject: [PATCH 09/58] remove #include Bidiagonalization, and add missing ";" --- Eigen/SVD | 1 - Eigen/src/Core/Assign.h | 30 +++++++++++++++--------------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/Eigen/SVD b/Eigen/SVD index 3aab35118..625071a75 100644 --- a/Eigen/SVD +++ b/Eigen/SVD @@ -22,7 +22,6 @@ namespace Eigen { * \endcode */ -#include "src/SVD/Bidiagonalization.h" #include "src/SVD/SVD.h" #include "src/SVD/JacobiSquareSVD.h" diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 68b6c817d..1d023286f 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -96,21 +96,21 @@ public: static void debug() { - EIGEN_DEBUG_VAR(DstIsAligned) - EIGEN_DEBUG_VAR(SrcIsAligned) - EIGEN_DEBUG_VAR(SrcAlignment) - EIGEN_DEBUG_VAR(InnerSize) - EIGEN_DEBUG_VAR(InnerMaxSize) - EIGEN_DEBUG_VAR(PacketSize) - EIGEN_DEBUG_VAR(MightVectorize) - EIGEN_DEBUG_VAR(MayInnerVectorize) - EIGEN_DEBUG_VAR(MayLinearVectorize) - EIGEN_DEBUG_VAR(MaySliceVectorize) - EIGEN_DEBUG_VAR(Vectorization) - EIGEN_DEBUG_VAR(UnrollingLimit) - EIGEN_DEBUG_VAR(MayUnrollCompletely) - EIGEN_DEBUG_VAR(MayUnrollInner) - EIGEN_DEBUG_VAR(Unrolling) + EIGEN_DEBUG_VAR(DstIsAligned); + EIGEN_DEBUG_VAR(SrcIsAligned); + EIGEN_DEBUG_VAR(SrcAlignment); + EIGEN_DEBUG_VAR(InnerSize); + EIGEN_DEBUG_VAR(InnerMaxSize); + EIGEN_DEBUG_VAR(PacketSize); + EIGEN_DEBUG_VAR(MightVectorize); + EIGEN_DEBUG_VAR(MayInnerVectorize); + EIGEN_DEBUG_VAR(MayLinearVectorize); + EIGEN_DEBUG_VAR(MaySliceVectorize); + EIGEN_DEBUG_VAR(Vectorization); + EIGEN_DEBUG_VAR(UnrollingLimit); + EIGEN_DEBUG_VAR(MayUnrollCompletely); + EIGEN_DEBUG_VAR(MayUnrollInner); + EIGEN_DEBUG_VAR(Unrolling); } }; From a4f664251863907604d43be70a41cc4c1dddd42a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 11 Aug 2009 15:11:47 +0200 Subject: [PATCH 10/58] fix issue #36 (missing return *this in Rotation2D --- Eigen/src/Geometry/Rotation2D.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Geometry/Rotation2D.h b/Eigen/src/Geometry/Rotation2D.h index 594d5272b..b80fcace0 100644 --- a/Eigen/src/Geometry/Rotation2D.h +++ b/Eigen/src/Geometry/Rotation2D.h @@ -85,7 +85,7 @@ public: /** Concatenates two rotations */ inline Rotation2D& operator*=(const Rotation2D& other) - { return m_angle += other.m_angle; } + { return m_angle += other.m_angle; return *this; } /** Applies the rotation to a 2D vector */ Vector2 operator* (const Vector2& vec) const From afbd73b5cdc1ce9b8ae54e9dd08332c870cf54d2 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 11 Aug 2009 15:15:06 +0200 Subject: [PATCH 11/58] overload operartor* with a ProductBase such that "scalar * (mat * mat)" is optimized as one could naturally expect --- Eigen/src/Core/Product.h | 6 +- Eigen/src/Core/ProductBase.h | 69 +++++++++++++++++-- Eigen/src/Core/products/GeneralMatrixMatrix.h | 2 +- .../Core/products/SelfadjointMatrixMatrix.h | 2 +- .../Core/products/SelfadjointMatrixVector.h | 2 +- .../Core/products/TriangularMatrixMatrix.h | 2 +- .../Core/products/TriangularMatrixVector.h | 2 +- test/product_notemporary.cpp | 2 +- test/product_symm.cpp | 5 ++ 9 files changed, 76 insertions(+), 16 deletions(-) diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 18f14f75a..610d5c84a 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -153,7 +153,7 @@ class GeneralProduct GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} - template void addTo(Dest& dst, Scalar alpha) const + template void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(dst.rows()==1 && dst.cols()==1); dst.coeffRef(0,0) += alpha * (m_lhs.cwise()*m_rhs).sum(); @@ -179,7 +179,7 @@ class GeneralProduct GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} - template void addTo(Dest& dest, Scalar alpha) const + template void scaleAndAddTo(Dest& dest, Scalar alpha) const { ei_outer_product_selector::run(*this, dest, alpha); } @@ -236,7 +236,7 @@ class GeneralProduct enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight }; typedef typename ei_meta_if::ret MatrixType; - template void addTo(Dest& dst, Scalar alpha) const + template void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols()); ei_gemv_selector inline int cols() const { return m_rhs.cols(); } template - inline void evalTo(Dest& dst) const { dst.setZero(); addTo(dst,1); } + inline void evalTo(Dest& dst) const { dst.setZero(); scaleAndAddTo(dst,1); } template - inline void addTo(Dest& dst) const { addTo(dst,1); } + inline void addTo(Dest& dst) const { scaleAndAddTo(dst,1); } template - inline void subTo(Dest& dst) const { addTo(dst,-1); } + inline void subTo(Dest& dst) const { scaleAndAddTo(dst,-1); } template - inline void addTo(Dest& dst,Scalar alpha) const { derived().addTo(dst,alpha); } + inline void scaleAndAddTo(Dest& dst,Scalar alpha) const { derived().scaleAndAddTo(dst,alpha); } PlainMatrixType eval() const { @@ -141,13 +141,68 @@ class ProductBase : public MatrixBase void coeffRef(int); }; +template +class ScaledProduct; + +// Note that these two operator* functions are not defined as member +// functions of ProductBase, because, otherwise we would have to +// define all overloads defined in MatrixBase. Furthermore, Using +// "using Base::operator*" would not work with MSVC. +template +const ScaledProduct operator*(const ProductBase& prod, typename Derived::Scalar x) +{ return ScaledProduct(prod.derived(), x); } + +template +const ScaledProduct operator*(typename Derived::Scalar x,const ProductBase& prod) +{ return ScaledProduct(prod.derived(), x); } + +template +struct ei_traits > + : ei_traits, + typename NestedProduct::_LhsNested, + typename NestedProduct::_RhsNested> > +{}; + +template +class ScaledProduct + : public ProductBase, + typename NestedProduct::_LhsNested, + typename NestedProduct::_RhsNested> +{ + public: + typedef ProductBase, + typename NestedProduct::_LhsNested, + typename NestedProduct::_RhsNested> Base; + typedef typename Base::Scalar Scalar; +// EIGEN_PRODUCT_PUBLIC_INTERFACE(ScaledProduct) + + ScaledProduct(const NestedProduct& prod, Scalar& x) + : Base(prod.lhs(),prod.rhs()), m_prod(prod), m_alpha(x) {} + + template + inline void evalTo(Dest& dst) const { dst.setZero(); scaleAndAddTo(dst,m_alpha); } + + template + inline void addTo(Dest& dst) const { scaleAndAddTo(dst,m_alpha); } + + template + inline void subTo(Dest& dst) const { scaleAndAddTo(dst,-m_alpha); } + + template + inline void scaleAndAddTo(Dest& dst,Scalar alpha) const { m_prod.derived().scaleAndAddTo(dst,alpha); } + + protected: + const NestedProduct& m_prod; + Scalar m_alpha; +}; + /** \internal * Overloaded to perform an efficient C = (A*B).lazy() */ template template Derived& MatrixBase::lazyAssign(const ProductBase& other) { - other.evalTo(derived()); return derived(); + other.derived().evalTo(derived()); return derived(); } /** \internal @@ -157,7 +212,7 @@ template Derived& MatrixBase::operator+=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) { - other._expression().addTo(derived()); return derived(); + other._expression().derived().addTo(derived()); return derived(); } /** \internal @@ -167,7 +222,7 @@ template Derived& MatrixBase::operator-=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) { - other._expression().subTo(derived()); return derived(); + other._expression().derived().subTo(derived()); return derived(); } #endif // EIGEN_PRODUCTBASE_H diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index ff0f2c1b4..8b3b13266 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -137,7 +137,7 @@ class GeneralProduct GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} - template void addTo(Dest& dst, Scalar alpha) const + template void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols()); diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h index 358da3752..5e025b90b 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h +++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h @@ -375,7 +375,7 @@ struct SelfadjointProductMatrix RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit }; - template void addTo(Dest& dst, Scalar alpha) const + template void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols()); diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index f0004cdb9..c2c33d5b8 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -175,7 +175,7 @@ struct SelfadjointProductMatrix SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} - template void addTo(Dest& dst, Scalar alpha) const + template void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols()); diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h index c2ee39e79..701ccb644 100644 --- a/Eigen/src/Core/products/TriangularMatrixMatrix.h +++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h @@ -333,7 +333,7 @@ struct TriangularProduct TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} - template void addTo(Dest& dst, Scalar alpha) const + template void scaleAndAddTo(Dest& dst, Scalar alpha) const { const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs); const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs); diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index a21afa2f6..620b090b9 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -130,7 +130,7 @@ struct TriangularProduct TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} - template void addTo(Dest& dst, Scalar alpha) const + template void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols()); diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index f4311e495..d5d996e49 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -79,7 +79,7 @@ template void product_notemporary(const MatrixType& m) VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * (m1*s3+m2*s2).adjoint()).lazy(), 1); VERIFY_EVALUATION_COUNT( m3 = ((s1 * m1).adjoint() * s2 * m2).lazy(), 0); VERIFY_EVALUATION_COUNT( m3 -= (s1 * (-m1*s3).adjoint() * (s2 * m2 * s3)).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 -= (s1 * (m1.transpose() * m2)).lazy(), 1); + VERIFY_EVALUATION_COUNT( m3 -= (s1 * (m1.transpose() * m2)).lazy(), 0); VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) += (-m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint()).lazy() ), 0); VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) -= (s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1)).lazy() ), 0); diff --git a/test/product_symm.cpp b/test/product_symm.cpp index 88bac878b..1300928a2 100644 --- a/test/product_symm.cpp +++ b/test/product_symm.cpp @@ -94,6 +94,11 @@ template void symm(int size = Size, in VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView() * (s2*rhs3).conjugate(), rhs13 = (s1*m1.adjoint()) * (s2*rhs3).conjugate()); + + m2 = m1.template triangularView(); rhs13 = rhs12; + VERIFY_IS_APPROX(rhs12 += (s1 * ((m2.adjoint()).template selfadjointView() * (s2*rhs3).conjugate())).lazy(), + rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate()); + // test matrix * selfadjoint symm_extra::run(m1,m2,rhs2,rhs22,rhs23,s1,s2); From ce033ebdfe817a21e648cad3f2bb6db76cb8fa6a Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Tue, 11 Aug 2009 16:12:34 -0400 Subject: [PATCH 12/58] add EIGEN_DEBUG_VAR --- Eigen/src/Core/Assign.h | 30 +++++++++++++++--------------- Eigen/src/Core/util/Macros.h | 2 ++ 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 1d023286f..4bd1046a7 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -96,21 +96,21 @@ public: static void debug() { - EIGEN_DEBUG_VAR(DstIsAligned); - EIGEN_DEBUG_VAR(SrcIsAligned); - EIGEN_DEBUG_VAR(SrcAlignment); - EIGEN_DEBUG_VAR(InnerSize); - EIGEN_DEBUG_VAR(InnerMaxSize); - EIGEN_DEBUG_VAR(PacketSize); - EIGEN_DEBUG_VAR(MightVectorize); - EIGEN_DEBUG_VAR(MayInnerVectorize); - EIGEN_DEBUG_VAR(MayLinearVectorize); - EIGEN_DEBUG_VAR(MaySliceVectorize); - EIGEN_DEBUG_VAR(Vectorization); - EIGEN_DEBUG_VAR(UnrollingLimit); - EIGEN_DEBUG_VAR(MayUnrollCompletely); - EIGEN_DEBUG_VAR(MayUnrollInner); - EIGEN_DEBUG_VAR(Unrolling); + EIGEN_DEBUG_VAR(DstIsAligned) + EIGEN_DEBUG_VAR(SrcIsAligned) + EIGEN_DEBUG_VAR(SrcAlignment) + EIGEN_DEBUG_VAR(InnerSize) + EIGEN_DEBUG_VAR(InnerMaxSize) + EIGEN_DEBUG_VAR(PacketSize) + EIGEN_DEBUG_VAR(MightVectorize) + EIGEN_DEBUG_VAR(MayInnerVectorize) + EIGEN_DEBUG_VAR(MayLinearVectorize) + EIGEN_DEBUG_VAR(MaySliceVectorize) + EIGEN_DEBUG_VAR(Vectorization) + EIGEN_DEBUG_VAR(UnrollingLimit) + EIGEN_DEBUG_VAR(MayUnrollCompletely) + EIGEN_DEBUG_VAR(MayUnrollInner) + EIGEN_DEBUG_VAR(Unrolling) } }; diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 956b609a0..f3f2b9066 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -111,6 +111,8 @@ #define EIGEN_FAST_MATH 1 #endif +#define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl; + #define USING_PART_OF_NAMESPACE_EIGEN \ EIGEN_USING_MATRIX_TYPEDEFS \ using Eigen::Matrix; \ From 22d65d47d085a9b693d7777f646a87cb3d51a06a Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Wed, 12 Aug 2009 02:35:07 -0400 Subject: [PATCH 13/58] finally, the good approach was two-sided Jacobi. Indeed, it allows to guarantee the precision of the output, which is very valuable. Here, we guarantee that the diagonal matrix returned by the SVD is actually diagonal, to machine precision. Performance isn't bad at all at 50% of the current householder SVD performance for a 200x200 matrix (no vectorization) and we have lots of room for improvement. --- Eigen/src/Core/MatrixBase.h | 5 ++- Eigen/src/Jacobi/CMakeLists.txt | 6 +++ Eigen/src/Jacobi/Jacobi.h | 30 ++++++++----- Eigen/src/SVD/JacobiSquareSVD.h | 74 ++++++++++++++++----------------- 4 files changed, 63 insertions(+), 52 deletions(-) create mode 100644 Eigen/src/Jacobi/CMakeLists.txt diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 688b7c4d0..52d5f680c 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -782,8 +782,9 @@ template class MatrixBase void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s); void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s); - bool makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s); - bool makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s); + bool makeJacobi(int p, int q, Scalar *c, Scalar *s) const; + bool makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const; + bool makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const; #ifdef EIGEN_MATRIXBASE_PLUGIN #include EIGEN_MATRIXBASE_PLUGIN diff --git a/Eigen/src/Jacobi/CMakeLists.txt b/Eigen/src/Jacobi/CMakeLists.txt new file mode 100644 index 000000000..490dac626 --- /dev/null +++ b/Eigen/src/Jacobi/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_Jacobi_SRCS "*.h") + +INSTALL(FILES + ${Eigen_Jacobi_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Jacobi COMPONENT Devel + ) diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 993a723ab..5866ac44f 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -48,13 +48,13 @@ void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s } template -bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar max_coeff, Scalar *c, Scalar *s) +bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) { - if(ei_abs(y) < max_coeff * 0.5 * machine_epsilon()) + if(ei_abs(y) < ei_abs(z-x) * 0.5 * machine_epsilon()) { *c = Scalar(1); *s = Scalar(0); - return true; + return false; } else { @@ -67,23 +67,31 @@ bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar max_coeff, Scalar *c, Sc t = Scalar(1) / (tau - w); *c = Scalar(1) / ei_sqrt(1 + ei_abs2(t)); *s = *c * t; - return false; + return true; } } template -inline bool MatrixBase::makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s) +inline bool MatrixBase::makeJacobi(int p, int q, Scalar *c, Scalar *s) const { - return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), max_coeff, c, s); + return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), c, s); } template -inline bool MatrixBase::makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s) +inline bool MatrixBase::makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const { - return ei_makeJacobi(col(p).squaredNorm(), - col(p).dot(col(q)), - col(q).squaredNorm(), - max_coeff, + return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(q,p)), + ei_conj(coeff(p,p))*coeff(p,q) + ei_conj(coeff(q,p))*coeff(q,q), + ei_abs2(coeff(p,q)) + ei_abs2(coeff(q,q)), + c,s); +} + +template +inline bool MatrixBase::makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const +{ + return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(p,q)), + ei_conj(coeff(q,p))*coeff(p,p) + ei_conj(coeff(q,q))*coeff(p,q), + ei_abs2(coeff(q,p)) + ei_abs2(coeff(q,q)), c,s); } diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h index ad55735fc..18416e312 100644 --- a/Eigen/src/SVD/JacobiSquareSVD.h +++ b/Eigen/src/SVD/JacobiSquareSVD.h @@ -102,69 +102,65 @@ void JacobiSquareSVD::compute(const MatrixType& if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); m_singularValues.resize(size); - RealScalar max_coeff = work_matrix.cwise().abs().maxCoeff(); - for(int k = 1; k < 40; ++k) { + while(true) + { bool finished = true; for(int p = 1; p < size; ++p) { for(int q = 0; q < p; ++q) { Scalar c, s; - finished &= work_matrix.makeJacobiForAtA(p,q,max_coeff,&c,&s); - work_matrix.applyJacobiOnTheRight(p,q,c,s); - if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s); + if(work_matrix.makeJacobiForAtA(p,q,&c,&s)) + { + work_matrix.applyJacobiOnTheRight(p,q,c,s); + if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s); + } + if(work_matrix.makeJacobiForAAt(p,q,&c,&s)) + { + work_matrix.applyJacobiOnTheLeft(p,q,c,s); + if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s); + if(std::max(ei_abs(work_matrix.coeff(p,q)), ei_abs(work_matrix.coeff(q,p))) > std::max(ei_abs(work_matrix.coeff(q,q)), ei_abs(work_matrix.coeff(p,p)) )) + { + work_matrix.row(p).swap(work_matrix.row(q)); + if(ComputeU) m_matrixU.col(p).swap(m_matrixU.col(q)); + } + } + } + } + RealScalar biggest = work_matrix.diagonal().cwise().abs().maxCoeff(); + for(int p = 0; p < size; ++p) + { + for(int q = 0; q < size; ++q) + { + if(p!=q && ei_abs(work_matrix.coeff(p,q)) > biggest * machine_epsilon()) finished = false; } } if(finished) break; } - + + m_singularValues = work_matrix.diagonal().cwise().abs(); + RealScalar biggestSingularValue = m_singularValues.maxCoeff(); + for(int i = 0; i < size; ++i) { - m_singularValues.coeffRef(i) = work_matrix.col(i).norm(); + RealScalar a = ei_abs(work_matrix.coeff(i,i)); + m_singularValues.coeffRef(i) = a; + if(ComputeU && !ei_isMuchSmallerThan(a, biggestSingularValue)) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a; } - int first_zero = size; - RealScalar biggest = m_singularValues.maxCoeff(); for(int i = 0; i < size; i++) { int pos; - RealScalar biggest_remaining = m_singularValues.end(size-i).maxCoeff(&pos); - if(first_zero == size && ei_isMuchSmallerThan(biggest_remaining, biggest)) first_zero = pos + i; + m_singularValues.end(size-i).maxCoeff(&pos); if(pos) { pos += i; std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); - if(ComputeU) work_matrix.col(pos).swap(work_matrix.col(i)); + if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); } } - - if(ComputeU) - { - for(int i = 0; i < first_zero; ++i) - { - m_matrixU.col(i) = work_matrix.col(i) / m_singularValues.coeff(i); - } - if(first_zero < size) - { - for(int i = first_zero; i < size; ++i) - { - for(int j = 0; j < size; ++j) - { - m_matrixU.col(i).setZero(); - m_matrixU.coeffRef(j,i) = Scalar(1); - for(int k = 0; k < first_zero; ++k) - m_matrixU.col(i) -= m_matrixU.col(i).dot(m_matrixU.col(k)) * m_matrixU.col(k); - RealScalar n = m_matrixU.col(i).norm(); - if(!ei_isMuchSmallerThan(n, biggest)) - { - m_matrixU.col(i) /= n; - break; - } - } - } - } - } + m_isInitialized = true; } #endif // EIGEN_JACOBISQUARESVD_H From 309d540d4a71527a4dba778dc5221b81fd18f540 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Wed, 12 Aug 2009 10:14:15 -0400 Subject: [PATCH 14/58] add parentheses; hopefully this solves Koldos MSVC compilation issue... --- Eigen/src/Core/BandMatrix.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h index c22696992..d3997397d 100644 --- a/Eigen/src/Core/BandMatrix.h +++ b/Eigen/src/Core/BandMatrix.h @@ -130,14 +130,14 @@ class BandMatrix : public AnyMatrixBase struct DiagonalIntReturnType { enum { - ReturnOpposite = (Options&SelfAdjoint) && (Index>0 && Supers==0 || Index<0 && Subs==0), + ReturnOpposite = (Options&SelfAdjoint) && ((Index>0 && Supers==0) || (Index<0 && Subs==0)), Conjugate = ReturnOpposite && NumTraits::IsComplex, ActualIndex = ReturnOpposite ? -Index : Index, - DiagonalSize = RowsAtCompileTime==Dynamic || ColsAtCompileTime==Dynamic + DiagonalSize = (RowsAtCompileTime==Dynamic || ColsAtCompileTime==Dynamic) ? Dynamic - : ActualIndex<0 + : (ActualIndex<0 ? EIGEN_ENUM_MIN(ColsAtCompileTime, RowsAtCompileTime + ActualIndex) - : EIGEN_ENUM_MIN(RowsAtCompileTime, ColsAtCompileTime - ActualIndex) + : EIGEN_ENUM_MIN(RowsAtCompileTime, ColsAtCompileTime - ActualIndex)) }; typedef Block BuildType; typedef typename ei_meta_if Date: Wed, 12 Aug 2009 15:44:22 +0100 Subject: [PATCH 15/58] Add support for matrix exponential of floats and complex numbers. --- .../src/MatrixFunctions/MatrixExponential.h | 331 +++++++++++++----- unsupported/test/matrixExponential.cpp | 64 +++- 2 files changed, 291 insertions(+), 104 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 39d2809d4..dd25d7f3d 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -25,13 +25,9 @@ #ifndef EIGEN_MATRIX_EXPONENTIAL #define EIGEN_MATRIX_EXPONENTIAL -#ifdef _MSC_VER -template Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); } -#endif - -/** Compute the matrix exponential. +/** \brief Compute the matrix exponential. * - * \param M matrix whose exponential is to be computed. + * \param M matrix whose exponential is to be computed. * \param result pointer to the matrix in which to store the result. * * The matrix exponential of \f$ M \f$ is defined by @@ -58,103 +54,264 @@ template Scalar log2(Scalar v) { return std::log(v)/std::log(S * SIAM J. %Matrix Anal. Applic., 26:1179–1193, * 2005. * - * \note Currently, \p M has to be a matrix of \c double . + * \note \p M has to be a matrix of \c float, \c double, + * \c complex or \c complex . */ template -void ei_matrix_exponential(const MatrixBase &M, typename ei_plain_matrix_type::type* result) -{ - typedef typename ei_traits::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef typename ei_plain_matrix_type::type PlainMatrixType; +EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase &M, + typename MatrixBase::PlainMatrixType* result); - ei_assert(M.rows() == M.cols()); - EIGEN_STATIC_ASSERT(NumTraits::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) - PlainMatrixType num, den, U, V; - PlainMatrixType Id = PlainMatrixType::Identity(M.rows(), M.cols()); - typename ei_eval::type Meval = M.eval(); - RealScalar l1norm = Meval.cwise().abs().colwise().sum().maxCoeff(); - int squarings = 0; - - // Choose degree of Pade approximant, depending on norm of M - if (l1norm < 1.495585217958292e-002) { - - // Use (3,3)-Pade +/** \internal \brief Internal helper functions for computing the + * matrix exponential. + */ +namespace MatrixExponentialInternal { + +#ifdef _MSC_VER + template Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); } +#endif + + /** \internal \brief Compute the (3,3)-Padé approximant to + * the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. + * + * \param M Argument of matrix exponential + * \param Id Identity matrix of same size as M + * \param tmp Temporary storage, to be provided by the caller + * \param M2 Temporary storage, to be provided by the caller + * \param U Even-degree terms in numerator of Padé approximant + * \param V Odd-degree terms in numerator of Padé approximant + */ + template + EIGEN_STRONG_INLINE void pade3(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, + MatrixType& M2, MatrixType& U, MatrixType& V) + { + typedef typename ei_traits::Scalar Scalar; const Scalar b[] = {120., 60., 12., 1.}; - PlainMatrixType M2; - M2 = (Meval * Meval).lazy(); - num = b[3]*M2 + b[1]*Id; - U = (Meval * num).lazy(); + M2 = (M * M).lazy(); + tmp = b[3]*M2 + b[1]*Id; + U = (M * tmp).lazy(); V = b[2]*M2 + b[0]*Id; - - } else if (l1norm < 2.539398330063230e-001) { - - // Use (5,5)-Pade + } + + /** \internal \brief Compute the (5,5)-Padé approximant to + * the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. + * + * \param M Argument of matrix exponential + * \param Id Identity matrix of same size as M + * \param tmp Temporary storage, to be provided by the caller + * \param M2 Temporary storage, to be provided by the caller + * \param U Even-degree terms in numerator of Padé approximant + * \param V Odd-degree terms in numerator of Padé approximant + */ + template + EIGEN_STRONG_INLINE void pade5(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, + MatrixType& M2, MatrixType& U, MatrixType& V) + { + typedef typename ei_traits::Scalar Scalar; const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.}; - PlainMatrixType M2, M4; - M2 = (Meval * Meval).lazy(); - M4 = (M2 * M2).lazy(); - num = b[5]*M4 + b[3]*M2 + b[1]*Id; - U = (Meval * num).lazy(); + M2 = (M * M).lazy(); + MatrixType M4 = (M2 * M2).lazy(); + tmp = b[5]*M4 + b[3]*M2 + b[1]*Id; + U = (M * tmp).lazy(); V = b[4]*M4 + b[2]*M2 + b[0]*Id; - - } else if (l1norm < 9.504178996162932e-001) { - - // Use (7,7)-Pade + } + + /** \internal \brief Compute the (7,7)-Padé approximant to + * the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. + * + * \param M Argument of matrix exponential + * \param Id Identity matrix of same size as M + * \param tmp Temporary storage, to be provided by the caller + * \param M2 Temporary storage, to be provided by the caller + * \param U Even-degree terms in numerator of Padé approximant + * \param V Odd-degree terms in numerator of Padé approximant + */ + template + EIGEN_STRONG_INLINE void pade7(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, + MatrixType& M2, MatrixType& U, MatrixType& V) + { + typedef typename ei_traits::Scalar Scalar; const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.}; - PlainMatrixType M2, M4, M6; - M2 = (Meval * Meval).lazy(); - M4 = (M2 * M2).lazy(); - M6 = (M4 * M2).lazy(); - num = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U = (Meval * num).lazy(); + M2 = (M * M).lazy(); + MatrixType M4 = (M2 * M2).lazy(); + MatrixType M6 = (M4 * M2).lazy(); + tmp = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; + U = (M * tmp).lazy(); V = b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; - - } else if (l1norm < 2.097847961257068e+000) { - - // Use (9,9)-Pade + } + + /** \internal \brief Compute the (9,9)-Padé approximant to + * the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. + * + * \param M Argument of matrix exponential + * \param Id Identity matrix of same size as M + * \param tmp Temporary storage, to be provided by the caller + * \param M2 Temporary storage, to be provided by the caller + * \param U Even-degree terms in numerator of Padé approximant + * \param V Odd-degree terms in numerator of Padé approximant + */ + template + EIGEN_STRONG_INLINE void pade9(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, + MatrixType& M2, MatrixType& U, MatrixType& V) + { + typedef typename ei_traits::Scalar Scalar; const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240., - 2162160., 110880., 3960., 90., 1.}; - PlainMatrixType M2, M4, M6, M8; - M2 = (Meval * Meval).lazy(); - M4 = (M2 * M2).lazy(); - M6 = (M4 * M2).lazy(); - M8 = (M6 * M2).lazy(); - num = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U = (Meval * num).lazy(); + 2162160., 110880., 3960., 90., 1.}; + M2 = (M * M).lazy(); + MatrixType M4 = (M2 * M2).lazy(); + MatrixType M6 = (M4 * M2).lazy(); + MatrixType M8 = (M6 * M2).lazy(); + tmp = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; + U = (M * tmp).lazy(); V = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; - - } else { - - // Use (13,13)-Pade; scale matrix by power of 2 so that its norm - // is small enough - - const Scalar maxnorm = 5.371920351148152; + } + + /** \internal \brief Compute the (13,13)-Padé approximant to + * the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. + * + * \param M Argument of matrix exponential + * \param Id Identity matrix of same size as M + * \param tmp Temporary storage, to be provided by the caller + * \param M2 Temporary storage, to be provided by the caller + * \param U Even-degree terms in numerator of Padé approximant + * \param V Odd-degree terms in numerator of Padé approximant + */ + template + EIGEN_STRONG_INLINE void pade13(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, + MatrixType& M2, MatrixType& U, MatrixType& V) + { + typedef typename ei_traits::Scalar Scalar; const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., - 1187353796428800., 129060195264000., 10559470521600., 670442572800., - 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.}; - - squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm))); - PlainMatrixType A, A2, A4, A6; - A = Meval / pow(Scalar(2), squarings); - A2 = (A * A).lazy(); - A4 = (A2 * A2).lazy(); - A6 = (A4 * A2).lazy(); - num = b[13]*A6 + b[11]*A4 + b[9]*A2; - V = (A6 * num).lazy(); - num = V + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*Id; - U = (A * num).lazy(); - num = b[12]*A6 + b[10]*A4 + b[8]*A2; - V = (A6 * num).lazy() + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*Id; + 1187353796428800., 129060195264000., 10559470521600., 670442572800., + 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.}; + M2 = (M * M).lazy(); + MatrixType M4 = (M2 * M2).lazy(); + MatrixType M6 = (M4 * M2).lazy(); + V = b[13]*M6 + b[11]*M4 + b[9]*M2; + tmp = (M6 * V).lazy(); + tmp += b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; + U = (M * tmp).lazy(); + tmp = b[12]*M6 + b[10]*M4 + b[8]*M2; + V = (M6 * tmp).lazy(); + V += b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; + } + + /** \internal \brief Helper class for computing Padé + * approximants to the exponential. + */ + template ::Scalar>::Real> + struct computeUV_selector + { + /** \internal \brief Compute Padé approximant to the exponential. + * + * Computes \p U, \p V and \p squarings such that \f$ (V+U)(V-U)^{-1} \f$ + * is a Padé of \f$ \exp(2^{-\mbox{squarings}}M) \f$ + * around \f$ M = 0 \f$. The degree of the Padé + * approximant and the value of squarings are chosen such that + * the approximation error is no more than the round-off error. + * + * \param M Argument of matrix exponential + * \param Id Identity matrix of same size as M + * \param tmp1 Temporary storage, to be provided by the caller + * \param tmp2 Temporary storage, to be provided by the caller + * \param U Even-degree terms in numerator of Padé approximant + * \param V Odd-degree terms in numerator of Padé approximant + * \param l1norm L1 norm of M + * \param squarings Pointer to integer containing number of times + * that the result needs to be squared to find the + * matrix exponential + */ + static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, + MatrixType& U, MatrixType& V, float l1norm, int* squarings); + }; + + template + struct computeUV_selector + { + static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, + MatrixType& U, MatrixType& V, float l1norm, int* squarings) + { + *squarings = 0; + if (l1norm < 4.258730016922831e-001) { + pade3(M, Id, tmp1, tmp2, U, V); + } else if (l1norm < 1.880152677804762e+000) { + pade5(M, Id, tmp1, tmp2, U, V); + } else { + const float maxnorm = 3.925724783138660; + *squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm))); + MatrixType A = M / std::pow(typename ei_traits::Scalar(2), *squarings); + pade7(A, Id, tmp1, tmp2, U, V); + } + } + }; + + template + struct computeUV_selector + { + static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, + MatrixType& U, MatrixType& V, float l1norm, int* squarings) + { + *squarings = 0; + if (l1norm < 1.495585217958292e-002) { + pade3(M, Id, tmp1, tmp2, U, V); + } else if (l1norm < 2.539398330063230e-001) { + pade5(M, Id, tmp1, tmp2, U, V); + } else if (l1norm < 9.504178996162932e-001) { + pade7(M, Id, tmp1, tmp2, U, V); + } else if (l1norm < 2.097847961257068e+000) { + pade9(M, Id, tmp1, tmp2, U, V); + } else { + const double maxnorm = 5.371920351148152; + *squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm))); + MatrixType A = M / std::pow(typename ei_traits::Scalar(2), *squarings); + pade13(A, Id, tmp1, tmp2, U, V); + } + } + }; + + /** \internal \brief Compute the matrix exponential. + * + * \param M matrix whose exponential is to be computed. + * \param result pointer to the matrix in which to store the result. + */ + template + void compute(const MatrixType &M, MatrixType* result) + { + MatrixType num, den, U, V; + MatrixType Id = MatrixType::Identity(M.rows(), M.cols()); + float l1norm = M.cwise().abs().colwise().sum().maxCoeff(); + int squarings; + computeUV_selector::run(M, Id, num, den, U, V, l1norm, &squarings); + num = U + V; // numerator of Pade approximant + den = -U + V; // denominator of Pade approximant + den.partialLu().solve(num, result); + for (int i=0; i +EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase &M, + typename MatrixBase::PlainMatrixType* result) +{ + ei_assert(M.rows() == M.cols()); + MatrixExponentialInternal::compute(M.eval(), result); } #endif // EIGEN_MATRIX_EXPONENTIAL diff --git a/unsupported/test/matrixExponential.cpp b/unsupported/test/matrixExponential.cpp index 0ec4fae90..9a6c335d8 100644 --- a/unsupported/test/matrixExponential.cpp +++ b/unsupported/test/matrixExponential.cpp @@ -34,26 +34,47 @@ double binom(int n, int k) return res; } -void test2dRotation() +template +void test2dRotation(double tol) { - Matrix2d A, B, C; - double angle; + Matrix A, B, C; + T angle; + A << 0, 1, -1, 0; for (int i=0; i<=20; i++) { angle = pow(10, i / 5. - 2); - A << 0, angle, -angle, 0; B << cos(angle), sin(angle), -sin(angle), cos(angle); - ei_matrix_exponential(A, &C); - VERIFY(C.isApprox(B, 1e-14)); + ei_matrix_exponential(angle*A, &C); + VERIFY(C.isApprox(B, tol)); } } -void testPascal() +template +void test2dHyperbolicRotation(double tol) +{ + Matrix,2,2> A, B, C; + std::complex imagUnit(0,1); + T angle, ch, sh; + + for (int i=0; i<=20; i++) + { + angle = (i-10) / 2.0; + ch = std::cosh(angle); + sh = std::sinh(angle); + A << 0, angle*imagUnit, -angle*imagUnit, 0; + B << ch, sh*imagUnit, -sh*imagUnit, ch; + ei_matrix_exponential(A, &C); + VERIFY(C.isApprox(B, tol)); + } +} + +template +void testPascal(double tol) { for (int size=1; size<20; size++) { - MatrixXd A(size,size), B(size,size), C(size,size); + Matrix A(size,size), B(size,size), C(size,size); A.setZero(); for (int i=0; i void randomTest(const MatrixType& m) +template +void randomTest(const MatrixType& m, double tol) { /* this test covers the following files: Inverse.h @@ -80,16 +102,24 @@ template void randomTest(const MatrixType& m) m1 = MatrixType::Random(rows, cols); ei_matrix_exponential(m1, &m2); ei_matrix_exponential(-m1, &m3); - VERIFY(identity.isApprox(m2 * m3, 1e-13)); + VERIFY(identity.isApprox(m2 * m3, tol)); } } void test_matrixExponential() { - CALL_SUBTEST(test2dRotation()); - CALL_SUBTEST(testPascal()); - CALL_SUBTEST(randomTest(Matrix2d())); - CALL_SUBTEST(randomTest(Matrix3d())); - CALL_SUBTEST(randomTest(Matrix4d())); - CALL_SUBTEST(randomTest(MatrixXd(8,8))); + CALL_SUBTEST(test2dRotation(1e-14)); + CALL_SUBTEST(test2dRotation(1e-5)); + CALL_SUBTEST(test2dHyperbolicRotation(1e-14)); + CALL_SUBTEST(test2dHyperbolicRotation(1e-5)); + CALL_SUBTEST(testPascal(1e-5)); + CALL_SUBTEST(testPascal(1e-14)); + CALL_SUBTEST(randomTest(Matrix2d(), 1e-13)); + CALL_SUBTEST(randomTest(Matrix(), 1e-13)); + CALL_SUBTEST(randomTest(Matrix4cd(), 1e-13)); + CALL_SUBTEST(randomTest(MatrixXd(8,8), 1e-13)); + CALL_SUBTEST(randomTest(Matrix2f(), 1e-4)); + CALL_SUBTEST(randomTest(Matrix3cf(), 1e-4)); + CALL_SUBTEST(randomTest(Matrix4f(), 1e-4)); + CALL_SUBTEST(randomTest(MatrixXf(8,8), 1e-4)); } From 2b618a2c1631548b00af3fe923e9db76aca71825 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Wed, 12 Aug 2009 18:23:39 -0400 Subject: [PATCH 16/58] make jacobi SVD more robust after experimenting with very nasty matrices... it turns out to be better to repeat the jacobi steps on a given (p,q) pair until it is diagonal to machine precision, before going to the next (p,q) pair. it's also an optimization as experiments show that in a majority of cases this allows to find out that the (p,q) pair is already diagonal to machine precision. --- Eigen/src/Jacobi/Jacobi.h | 2 +- Eigen/src/SVD/JacobiSquareSVD.h | 41 ++++++++++++++++++++------------- 2 files changed, 26 insertions(+), 17 deletions(-) diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 5866ac44f..48195c1ec 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -50,7 +50,7 @@ void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s template bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) { - if(ei_abs(y) < ei_abs(z-x) * 0.5 * machine_epsilon()) + if(ei_abs(y) <= ei_abs(z-x) * 0.5 * machine_epsilon()) { *c = Scalar(1); *s = Scalar(0); diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h index 18416e312..18c3db980 100644 --- a/Eigen/src/SVD/JacobiSquareSVD.h +++ b/Eigen/src/SVD/JacobiSquareSVD.h @@ -102,14 +102,16 @@ void JacobiSquareSVD::compute(const MatrixType& if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); m_singularValues.resize(size); - while(true) + +sweep_again: + for(int p = 1; p < size; ++p) { - bool finished = true; - for(int p = 1; p < size; ++p) + for(int q = 0; q < p; ++q) { - for(int q = 0; q < p; ++q) - { - Scalar c, s; + Scalar c, s; + while(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) + > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*machine_epsilon()) + { if(work_matrix.makeJacobiForAtA(p,q,&c,&s)) { work_matrix.applyJacobiOnTheRight(p,q,c,s); @@ -117,9 +119,13 @@ void JacobiSquareSVD::compute(const MatrixType& } if(work_matrix.makeJacobiForAAt(p,q,&c,&s)) { + Scalar x = ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(p,q)); + Scalar y = ei_conj(work_matrix.coeff(q,p))*work_matrix.coeff(p,p) + ei_conj(work_matrix.coeff(q,q))*work_matrix.coeff(p,q); + Scalar z = ei_abs2(work_matrix.coeff(q,p)) + ei_abs2(work_matrix.coeff(q,q)); work_matrix.applyJacobiOnTheLeft(p,q,c,s); if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s); - if(std::max(ei_abs(work_matrix.coeff(p,q)), ei_abs(work_matrix.coeff(q,p))) > std::max(ei_abs(work_matrix.coeff(q,q)), ei_abs(work_matrix.coeff(p,p)) )) + if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) + > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q))) ) { work_matrix.row(p).swap(work_matrix.row(q)); if(ComputeU) m_matrixU.col(p).swap(m_matrixU.col(q)); @@ -127,15 +133,18 @@ void JacobiSquareSVD::compute(const MatrixType& } } } - RealScalar biggest = work_matrix.diagonal().cwise().abs().maxCoeff(); - for(int p = 0; p < size; ++p) - { - for(int q = 0; q < size; ++q) - { - if(p!=q && ei_abs(work_matrix.coeff(p,q)) > biggest * machine_epsilon()) finished = false; - } - } - if(finished) break; + } + + RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff(); + RealScalar maxAllowedOffDiag = biggestOnDiag * machine_epsilon(); + for(int p = 0; p < size; ++p) + { + for(int q = 0; q < p; ++q) + if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) + goto sweep_again; + for(int q = p+1; q < size; ++q) + if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) + goto sweep_again; } m_singularValues = work_matrix.diagonal().cwise().abs(); From 99802094e471bb64cc3c3e3d0b2391ba2a1e75a2 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Wed, 12 Aug 2009 18:30:37 -0400 Subject: [PATCH 17/58] do without an empirical homemade formula that i wasn't comfortable about... turns out it's not needed anymore and removing it seems to only increase the precision --- Eigen/src/Jacobi/Jacobi.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 48195c1ec..785dd26ab 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -50,7 +50,7 @@ void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s template bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) { - if(ei_abs(y) <= ei_abs(z-x) * 0.5 * machine_epsilon()) + if(y == 0) { *c = Scalar(1); *s = Scalar(0); From 1d80f561ad48a96e9852dd8bfdbc400574aa67e5 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 13 Aug 2009 22:50:55 -0400 Subject: [PATCH 18/58] apply change discussed on the list : * new default precision "-1" means use the current stream precision * otherwise, save and restore the stream precision --- Eigen/src/Core/IO.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h index 38ac3eef3..7cd2b683a 100644 --- a/Eigen/src/Core/IO.h +++ b/Eigen/src/Core/IO.h @@ -33,7 +33,7 @@ enum { Raw, AlignCols }; * \brief Stores a set of parameters controlling the way matrices are printed * * List of available parameters: - * - \b precision number of digits for floating point values + * - \b precision number of digits for floating point values. The default value -1 means that the current stream precision is used. * - \b flags can be either Raw (default) or AlignCols which aligns all the columns * - \b coeffSeparator string printed between two coefficients of the same row * - \b rowSeparator string printed between two rows @@ -50,7 +50,7 @@ enum { Raw, AlignCols }; struct IOFormat { /** Default contructor, see class IOFormat for the meaning of the parameters */ - IOFormat(int _precision=4, int _flags=Raw, + IOFormat(int _precision=-1, int _flags=Raw, const std::string& _coeffSeparator = " ", const std::string& _rowSeparator = "\n", const std::string& _rowPrefix="", const std::string& _rowSuffix="", const std::string& _matPrefix="", const std::string& _matSuffix="") @@ -134,12 +134,12 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm for(int i = 0; i < m.rows(); ++i) { std::stringstream sstr; - sstr.precision(fmt.precision); + if(fmt.precision != -1) sstr.precision(fmt.precision); sstr << m.coeff(i,j); width = std::max(width, int(sstr.str().length())); } } - s.precision(fmt.precision); + std::streamsize old_precision = s.precision(fmt.precision); s << fmt.matPrefix; for(int i = 0; i < m.rows(); ++i) { @@ -159,6 +159,7 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm s << fmt.rowSeparator; } s << fmt.matSuffix; + s.precision(old_precision); return s; } From 1b257a7620d122a05f71b6ebd52e14fe5f1b18ef Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 13 Aug 2009 11:42:02 +0200 Subject: [PATCH 19/58] add an optimized "apply in place a rotation in the plane", and make Jacobi and SelfAdjointEigenSolver use it => ~ x1.75 speedup for JacobiSVD and x2 for SelfAdjointEigenSolver --- Eigen/Core | 1 + Eigen/src/Core/products/RotationInThePlane.h | 127 +++++++++++++++++++ Eigen/src/Jacobi/Jacobi.h | 18 +-- Eigen/src/QR/SelfAdjointEigenSolver.h | 17 +-- 4 files changed, 136 insertions(+), 27 deletions(-) create mode 100644 Eigen/src/Core/products/RotationInThePlane.h diff --git a/Eigen/Core b/Eigen/Core index 28ed09035..0a1477bbf 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -194,6 +194,7 @@ namespace Eigen { #include "src/Core/products/TriangularMatrixVector.h" #include "src/Core/products/TriangularMatrixMatrix.h" #include "src/Core/products/TriangularSolverMatrix.h" +#include "src/Core/products/RotationInThePlane.h" #include "src/Core/BandMatrix.h" } // namespace Eigen diff --git a/Eigen/src/Core/products/RotationInThePlane.h b/Eigen/src/Core/products/RotationInThePlane.h new file mode 100644 index 000000000..aa0b4d640 --- /dev/null +++ b/Eigen/src/Core/products/RotationInThePlane.h @@ -0,0 +1,127 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_ROTATION_IN_THE_PLANE_H +#define EIGEN_ROTATION_IN_THE_PLANE_H + +/********************************************************************** +* This file implement ... +**********************************************************************/ + +template +struct ei_apply_rotation_in_the_plane_selector; + +template +void ei_apply_rotation_in_the_plane(VectorX& x, VectorY& y, typename VectorX::Scalar c, typename VectorY::Scalar s) +{ + ei_assert(x.size() == y.size()); + int size = x.size(); + int incrx = size ==1 ? 1 : &x.coeffRef(1) - &x.coeffRef(0); + int incry = size ==1 ? 1 : &y.coeffRef(1) - &y.coeffRef(0); + if (incrx==1 && incry==1) + ei_apply_rotation_in_the_plane_selector + ::run(&x.coeffRef(0), &y.coeffRef(0), x.size(), c, s, 1, 1); + else + ei_apply_rotation_in_the_plane_selector + ::run(&x.coeffRef(0), &y.coeffRef(0), x.size(), c, s, incrx, incry); +} + +template +struct ei_apply_rotation_in_the_plane_selector +{ + static EIGEN_DONT_INLINE void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int incrx, int incry) + { + for(int i=0; i vectorization +template +struct ei_apply_rotation_in_the_plane_selector +{ + static EIGEN_DONT_INLINE void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int, int) + { + typedef typename ei_packet_traits::type Packet; + enum { PacketSize = ei_packet_traits::size, Peeling = 2 }; + int alignedStart = ei_alignmentOffset(y, size); + int alignedEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize); + + const Packet pc = ei_pset1(c); + const Packet ps = ei_pset1(s); + + for(int i=0; i void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) { - for(int i = 0; i < cols(); ++i) - { - Scalar tmp = coeff(p,i); - coeffRef(p,i) = c * tmp - s * coeff(q,i); - coeffRef(q,i) = s * tmp + c * coeff(q,i); - } + RowXpr x(row(p)); + RowXpr y(row(q)); + ei_apply_rotation_in_the_plane(x, y, c, s); } template void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) { - for(int i = 0; i < rows(); ++i) - { - Scalar tmp = coeff(i,p); - coeffRef(i,p) = c * tmp - s * coeff(i,q); - coeffRef(i,q) = s * tmp + c * coeff(i,q); - } + ColXpr x(col(p)); + ColXpr y(col(q)); + ei_apply_rotation_in_the_plane(x, y, c, s); } template diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index b003fb4d7..a8e89ba1a 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -378,23 +378,10 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st if (matrixQ) { #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR + ei_apply_rotation_in_the_plane_selector::run(matrixQ+k, matrixQ+k+1, n, c, s, n, n); #else - int kn = k*n; - int kn1 = (k+1)*n; + ei_apply_rotation_in_the_plane_selector::run(matrixQ+k*n, matrixQ+k*n+n, n, c, s, 1, 1); #endif - // let's do the product manually to avoid the need of temporaries... - for (int i=0; i Date: Thu, 13 Aug 2009 09:53:47 -0400 Subject: [PATCH 20/58] apply Koldo's workaround for MSVC bug --- Eigen/src/Core/BandMatrix.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h index d3997397d..841d1786a 100644 --- a/Eigen/src/Core/BandMatrix.h +++ b/Eigen/src/Core/BandMatrix.h @@ -130,7 +130,7 @@ class BandMatrix : public AnyMatrixBase struct DiagonalIntReturnType { enum { - ReturnOpposite = (Options&SelfAdjoint) && ((Index>0 && Supers==0) || (Index<0 && Subs==0)), + ReturnOpposite = (Options&SelfAdjoint) && (((Index)>0 && Supers==0) || ((Index)<0 && Subs==0)), Conjugate = ReturnOpposite && NumTraits::IsComplex, ActualIndex = ReturnOpposite ? -Index : Index, DiagonalSize = (RowsAtCompileTime==Dynamic || ColsAtCompileTime==Dynamic) From 76a3089a43752859893422bdd3d41874c465a44e Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 13 Aug 2009 09:56:53 -0400 Subject: [PATCH 21/58] oops, don't set the precision to -1 !! --- Eigen/src/Core/IO.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h index 7cd2b683a..53e02412f 100644 --- a/Eigen/src/Core/IO.h +++ b/Eigen/src/Core/IO.h @@ -139,7 +139,8 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm width = std::max(width, int(sstr.str().length())); } } - std::streamsize old_precision = s.precision(fmt.precision); + std::streamsize old_precision; + if(fmt.precision != -1) old_precision = s.precision(fmt.precision); s << fmt.matPrefix; for(int i = 0; i < m.rows(); ++i) { @@ -159,7 +160,7 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm s << fmt.rowSeparator; } s << fmt.matSuffix; - s.precision(old_precision); + if(fmt.precision != -1) s.precision(old_precision); return s; } From f2536416da990f12e98d01806331ad8d78545863 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 13 Aug 2009 14:56:39 -0400 Subject: [PATCH 22/58] * remove EIGEN_DONT_INLINE that harm performance for small sizes * normalize left Jacobi rotations to avoid having to swap rows * set precision to 2*machine_epsilon instead of machine_epsilon, we lose 1 bit of precision but gain between 10% and 100% speed, plus reduce the risk that some day we hit a bad matrix where it's impossible to approach machine precision --- Eigen/src/Core/products/RotationInThePlane.h | 4 ++-- Eigen/src/Jacobi/Jacobi.h | 16 ++++++++++++++-- Eigen/src/SVD/JacobiSquareSVD.h | 15 ++++----------- 3 files changed, 20 insertions(+), 15 deletions(-) diff --git a/Eigen/src/Core/products/RotationInThePlane.h b/Eigen/src/Core/products/RotationInThePlane.h index aa0b4d640..7fc62c675 100644 --- a/Eigen/src/Core/products/RotationInThePlane.h +++ b/Eigen/src/Core/products/RotationInThePlane.h @@ -50,7 +50,7 @@ void ei_apply_rotation_in_the_plane(VectorX& x, VectorY& y, typename VectorX::Sc template struct ei_apply_rotation_in_the_plane_selector { - static EIGEN_DONT_INLINE void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int incrx, int incry) + static void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int incrx, int incry) { for(int i=0; i template struct ei_apply_rotation_in_the_plane_selector { - static EIGEN_DONT_INLINE void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int, int) + static void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int, int) { typedef typename ei_packet_traits::type Packet; enum { PacketSize = ei_packet_traits::size, Peeling = 2 }; diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 10ae9eb8f..40181cd08 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -26,7 +26,7 @@ #define EIGEN_JACOBI_H template -void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) +inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) { RowXpr x(row(p)); RowXpr y(row(q)); @@ -34,7 +34,7 @@ void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) } template -void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) +inline void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) { ColXpr x(col(p)); ColXpr y(col(q)); @@ -89,5 +89,17 @@ inline bool MatrixBase::makeJacobiForAAt(int p, int q, Scalar *c, Scala c,s); } +template +inline void ei_normalizeJacobi(Scalar *c, Scalar *s, const Scalar& x, const Scalar& y) +{ + Scalar a = x * *c - y * *s; + Scalar b = x * *s + y * *c; + if(ei_abs(b)>ei_abs(a)) { + Scalar x = *c; + *c = -*s; + *s = x; + } +} + #endif // EIGEN_JACOBI_H diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h index 18c3db980..82133f7be 100644 --- a/Eigen/src/SVD/JacobiSquareSVD.h +++ b/Eigen/src/SVD/JacobiSquareSVD.h @@ -102,6 +102,7 @@ void JacobiSquareSVD::compute(const MatrixType& if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); m_singularValues.resize(size); + const RealScalar precision = 2 * machine_epsilon(); sweep_again: for(int p = 1; p < size; ++p) @@ -110,7 +111,7 @@ sweep_again: { Scalar c, s; while(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) - > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*machine_epsilon()) + > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) { if(work_matrix.makeJacobiForAtA(p,q,&c,&s)) { @@ -119,24 +120,16 @@ sweep_again: } if(work_matrix.makeJacobiForAAt(p,q,&c,&s)) { - Scalar x = ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(p,q)); - Scalar y = ei_conj(work_matrix.coeff(q,p))*work_matrix.coeff(p,p) + ei_conj(work_matrix.coeff(q,q))*work_matrix.coeff(p,q); - Scalar z = ei_abs2(work_matrix.coeff(q,p)) + ei_abs2(work_matrix.coeff(q,q)); + ei_normalizeJacobi(&c, &s, work_matrix.coeff(p,p), work_matrix.coeff(q,p)), work_matrix.applyJacobiOnTheLeft(p,q,c,s); if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s); - if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) - > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q))) ) - { - work_matrix.row(p).swap(work_matrix.row(q)); - if(ComputeU) m_matrixU.col(p).swap(m_matrixU.col(q)); - } } } } } RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff(); - RealScalar maxAllowedOffDiag = biggestOnDiag * machine_epsilon(); + RealScalar maxAllowedOffDiag = biggestOnDiag * precision; for(int p = 0; p < size; ++p) { for(int q = 0; q < p; ++q) From 13e95f7f68f67ee0b66cfa37d9cdc9fdd71eeca8 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 14 Aug 2009 00:17:14 +0200 Subject: [PATCH 23/58] optimize "apply Jacobi" for small sizes, and move it to Jacobi.h --- Eigen/QR | 1 + Eigen/src/Core/products/RotationInThePlane.h | 127 ------------------- Eigen/src/Jacobi/Jacobi.h | 97 ++++++++++++++ 3 files changed, 98 insertions(+), 127 deletions(-) delete mode 100644 Eigen/src/Core/products/RotationInThePlane.h diff --git a/Eigen/QR b/Eigen/QR index 5f36d0987..ecebd32e4 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -6,6 +6,7 @@ #include "src/Core/util/DisableMSVCWarnings.h" #include "Cholesky" +#include "Jacobi" // Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module #if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) diff --git a/Eigen/src/Core/products/RotationInThePlane.h b/Eigen/src/Core/products/RotationInThePlane.h deleted file mode 100644 index 7fc62c675..000000000 --- a/Eigen/src/Core/products/RotationInThePlane.h +++ /dev/null @@ -1,127 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Gael Guennebaud -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_ROTATION_IN_THE_PLANE_H -#define EIGEN_ROTATION_IN_THE_PLANE_H - -/********************************************************************** -* This file implement ... -**********************************************************************/ - -template -struct ei_apply_rotation_in_the_plane_selector; - -template -void ei_apply_rotation_in_the_plane(VectorX& x, VectorY& y, typename VectorX::Scalar c, typename VectorY::Scalar s) -{ - ei_assert(x.size() == y.size()); - int size = x.size(); - int incrx = size ==1 ? 1 : &x.coeffRef(1) - &x.coeffRef(0); - int incry = size ==1 ? 1 : &y.coeffRef(1) - &y.coeffRef(0); - if (incrx==1 && incry==1) - ei_apply_rotation_in_the_plane_selector - ::run(&x.coeffRef(0), &y.coeffRef(0), x.size(), c, s, 1, 1); - else - ei_apply_rotation_in_the_plane_selector - ::run(&x.coeffRef(0), &y.coeffRef(0), x.size(), c, s, incrx, incry); -} - -template -struct ei_apply_rotation_in_the_plane_selector -{ - static void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int incrx, int incry) - { - for(int i=0; i vectorization -template -struct ei_apply_rotation_in_the_plane_selector -{ - static void run(Scalar* x, Scalar* y, int size, Scalar c, Scalar s, int, int) - { - typedef typename ei_packet_traits::type Packet; - enum { PacketSize = ei_packet_traits::size, Peeling = 2 }; - int alignedStart = ei_alignmentOffset(y, size); - int alignedEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize); - - const Packet pc = ei_pset1(c); - const Packet ps = ei_pset1(s); - - for(int i=0; i +// Copyright (C) 2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -25,6 +26,9 @@ #ifndef EIGEN_JACOBI_H #define EIGEN_JACOBI_H +template +void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s); + template inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) { @@ -101,5 +105,98 @@ inline void ei_normalizeJacobi(Scalar *c, Scalar *s, const Scalar& x, const Scal } } +template +void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s) +{ + typedef typename VectorX::Scalar Scalar; + ei_assert(_x.size() == _y.size()); + int size = _x.size(); + int incrx = size ==1 ? 1 : &_x.coeffRef(1) - &_x.coeffRef(0); + int incry = size ==1 ? 1 : &_y.coeffRef(1) - &_y.coeffRef(0); + + Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0); + Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0); + + if (incrx==1 && incry==1) + { + // both vectors are sequentially stored in memory => vectorization + typedef typename ei_packet_traits::type Packet; + enum { PacketSize = ei_packet_traits::size, Peeling = 2 }; + + int alignedStart = ei_alignmentOffset(y, size); + int alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize; + + const Packet pc = ei_pset1(c); + const Packet ps = ei_pset1(s); + + for(int i=0; i Date: Fri, 14 Aug 2009 09:49:33 +0200 Subject: [PATCH 24/58] oops forgot to remove the #include in Core --- Eigen/Core | 1 - 1 file changed, 1 deletion(-) diff --git a/Eigen/Core b/Eigen/Core index 0a1477bbf..28ed09035 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -194,7 +194,6 @@ namespace Eigen { #include "src/Core/products/TriangularMatrixVector.h" #include "src/Core/products/TriangularMatrixMatrix.h" #include "src/Core/products/TriangularSolverMatrix.h" -#include "src/Core/products/RotationInThePlane.h" #include "src/Core/BandMatrix.h" } // namespace Eigen From 6373c3cd0014ad8e0a6952915c3c6342b57c1863 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 14 Aug 2009 13:49:29 +0200 Subject: [PATCH 25/58] oops bis, I forgot that SelfAdjointEigneSolver directly called the selector... --- Eigen/src/QR/SelfAdjointEigenSolver.h | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index a8e89ba1a..d07d28424 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -377,11 +377,8 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st // G only modifies the two columns k and k+1 if (matrixQ) { - #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR - ei_apply_rotation_in_the_plane_selector::run(matrixQ+k, matrixQ+k+1, n, c, s, n, n); - #else - ei_apply_rotation_in_the_plane_selector::run(matrixQ+k*n, matrixQ+k*n+n, n, c, s, 1, 1); - #endif + Map > q(matrixQ,n,n); + q.applyJacobiOnTheRight(k,k+1,c,s); } } } From 22ae236d4ea404122ff68966f1572a236f239335 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 15:12:32 -0400 Subject: [PATCH 26/58] machine_epsilon -> epsilon as wrapper around numeric_traits --- Eigen/src/Cholesky/LDLT.h | 2 +- Eigen/src/Core/MathFunctions.h | 15 +++++++-------- Eigen/src/LU/LU.h | 2 +- Eigen/src/SVD/JacobiSquareSVD.h | 2 +- 4 files changed, 10 insertions(+), 11 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index e652fd213..1cb1294ac 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -180,7 +180,7 @@ void LDLT::compute(const MatrixType& a) // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical // Algorithms" page 217, also by Higham. - cutoff = ei_abs(machine_epsilon() * size * biggest_in_corner); + cutoff = ei_abs(epsilon() * size * biggest_in_corner); m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; } diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index d68c483b2..16091caf0 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2006-2008 Benoit Jacob +// Copyright (C) 2006-2009 Benoit Jacob // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -25,8 +25,13 @@ #ifndef EIGEN_MATHFUNCTIONS_H #define EIGEN_MATHFUNCTIONS_H +template typename NumTraits::Real epsilon() +{ + return std::numeric_limits::Real>::epsilon(); +} + template inline typename NumTraits::Real precision(); -template inline typename NumTraits::Real machine_epsilon(); + template inline T ei_random(T a, T b); template inline T ei_random(); template inline T ei_random_amplitude() @@ -51,7 +56,6 @@ template inline typename NumTraits::Real ei_hypot(T x, T y) **************/ template<> inline int precision() { return 0; } -template<> inline int machine_epsilon() { return 0; } inline int ei_real(int x) { return x; } inline int& ei_real_ref(int& x) { return x; } inline int ei_imag(int) { return 0; } @@ -106,7 +110,6 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision()) **************/ template<> inline float precision() { return 1e-5f; } -template<> inline float machine_epsilon() { return 1.192e-07f; } inline float ei_real(float x) { return x; } inline float& ei_real_ref(float& x) { return x; } inline float ei_imag(float) { return 0.f; } @@ -154,7 +157,6 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision inline double precision() { return 1e-12; } -template<> inline double machine_epsilon() { return 2.220e-16; } inline double ei_real(double x) { return x; } inline double& ei_real_ref(double& x) { return x; } @@ -203,7 +205,6 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision inline float precision >() { return precision(); } -template<> inline float machine_epsilon >() { return machine_epsilon(); } inline float ei_real(const std::complex& x) { return std::real(x); } inline float ei_imag(const std::complex& x) { return std::imag(x); } inline float& ei_real_ref(std::complex& x) { return reinterpret_cast(&x)[0]; } @@ -240,7 +241,6 @@ inline bool ei_isApprox(const std::complex& a, const std::complex& **********************/ template<> inline double precision >() { return precision(); } -template<> inline double machine_epsilon >() { return machine_epsilon(); } inline double ei_real(const std::complex& x) { return std::real(x); } inline double ei_imag(const std::complex& x) { return std::imag(x); } inline double& ei_real_ref(std::complex& x) { return reinterpret_cast(&x)[0]; } @@ -278,7 +278,6 @@ inline bool ei_isApprox(const std::complex& a, const std::complex inline long double precision() { return precision(); } -template<> inline long double machine_epsilon() { return 1.084e-19l; } inline long double ei_real(long double x) { return x; } inline long double& ei_real_ref(long double& x) { return x; } inline long double ei_imag(long double) { return 0.; } diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index f36951fda..733fa0cbc 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -390,7 +390,7 @@ void LU::compute(const MatrixType& matrix) // 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. - m_precision = machine_epsilon() * size; + m_precision = epsilon() * size; IntColVectorType rows_transpositions(matrix.rows()); IntRowVectorType cols_transpositions(matrix.cols()); diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h index 82133f7be..32adb1301 100644 --- a/Eigen/src/SVD/JacobiSquareSVD.h +++ b/Eigen/src/SVD/JacobiSquareSVD.h @@ -102,7 +102,7 @@ void JacobiSquareSVD::compute(const MatrixType& if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); m_singularValues.resize(size); - const RealScalar precision = 2 * machine_epsilon(); + const RealScalar precision = 2 * epsilon(); sweep_again: for(int p = 1; p < size; ++p) From 16abc0ba7f90293289e0fc7866f07d185d500ba6 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 15:49:14 -0400 Subject: [PATCH 27/58] try to support 16 bit platforms... optimistic, but can't hurt --- Eigen/src/Core/util/Constants.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 98b3bd2e4..21e7f62c3 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2008-2009 Gael Guennebaud -// Copyright (C) 2006-2008 Benoit Jacob +// Copyright (C) 2007-2009 Benoit Jacob // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -39,10 +39,9 @@ * Also, disabling compiler warnings for integer overflow, sounds like a bad idea. * - It should be a prime number, because for example the old value 10000 led to bugs with 100x100 matrices. * - * If you wish to port Eigen to a platform where sizeof(int)==2, it is perfectly possible to set Dynamic to, say, 97. - * However, changing the value of Dynamic breaks the ABI, as Dynamic is often used as a template parameter for Matrix. + * Changing the value of Dynamic breaks the ABI, as Dynamic is often used as a template parameter for Matrix. */ -const int Dynamic = 33331; +const int Dynamic = sizeof(int) >= 4 ? 33331 : 101; /** This value means +Infinity; it is currently used only as the p parameter to MatrixBase::lpNorm(). * The value Infinity there means the L-infinity norm. From 2f74801ca4b1cd82cd95cd98ab1ee0857efdcac2 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 16:31:42 -0400 Subject: [PATCH 28/58] as discussed on list: default to align cols, reorganize parameters accordingly so that the default corresponds to 0 flag, and implement FullPrecision output (non-default). --- Eigen/src/Core/IO.h | 44 +++++++++++++++++++++++++++------- Eigen/src/Core/MathFunctions.h | 2 +- 2 files changed, 36 insertions(+), 10 deletions(-) diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h index 53e02412f..03e1af17d 100644 --- a/Eigen/src/Core/IO.h +++ b/Eigen/src/Core/IO.h @@ -26,15 +26,22 @@ #ifndef EIGEN_IO_H #define EIGEN_IO_H -enum { Raw, AlignCols }; +enum { DontAlignCols = 1 }; +enum { StreamPrecision = -1, + FullPrecision = -2 }; /** \class IOFormat * * \brief Stores a set of parameters controlling the way matrices are printed * * List of available parameters: - * - \b precision number of digits for floating point values. The default value -1 means that the current stream precision is used. - * - \b flags can be either Raw (default) or AlignCols which aligns all the columns + * - \b precision number of digits for floating point values, or one of the special constants \c StreamPrecision and \c FullPrecision. + * The default is the special value \c StreamPrecision which means to use the + * stream's own precision setting, as set for instance using \c cout.precision(3). The other special value + * \c FullPrecision means that the number of digits will be computed to match the full precision of each floating-point + * type. + * - \b flags an OR-ed combination of flags, the default value is 0, the only currently available flag is \c DontAlignCols which + * allows to disable the alignment of columns, resulting in faster code. * - \b coeffSeparator string printed between two coefficients of the same row * - \b rowSeparator string printed between two rows * - \b rowPrefix string printed at the beginning of each row @@ -50,7 +57,7 @@ enum { Raw, AlignCols }; struct IOFormat { /** Default contructor, see class IOFormat for the meaning of the parameters */ - IOFormat(int _precision=-1, int _flags=Raw, + IOFormat(int _precision = StreamPrecision, int _flags = 0, const std::string& _coeffSeparator = " ", const std::string& _rowSeparator = "\n", const std::string& _rowPrefix="", const std::string& _rowSuffix="", const std::string& _matPrefix="", const std::string& _matSuffix="") @@ -125,22 +132,41 @@ template std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOFormat& fmt) { const typename Derived::Nested m = _m; - + typedef typename Derived::Scalar Scalar; + int width = 0; - if (fmt.flags & AlignCols) + + std::streamsize explicit_precision; + if(fmt.precision == StreamPrecision) + { + explicit_precision = 0; + } + else if(fmt.precision == FullPrecision) + { + explicit_precision = NumTraits::HasFloatingPoint + ? std::ceil(-ei_log(epsilon())/ei_log(10.0)) + : 0; + } + else + { + explicit_precision = fmt.precision; + } + + bool align_cols = !(fmt.flags & DontAlignCols); + if(align_cols) { // compute the largest width for(int j = 1; j < m.cols(); ++j) for(int i = 0; i < m.rows(); ++i) { std::stringstream sstr; - if(fmt.precision != -1) sstr.precision(fmt.precision); + if(explicit_precision) sstr.precision(explicit_precision); sstr << m.coeff(i,j); width = std::max(width, int(sstr.str().length())); } } std::streamsize old_precision; - if(fmt.precision != -1) old_precision = s.precision(fmt.precision); + if(explicit_precision) old_precision = s.precision(explicit_precision); s << fmt.matPrefix; for(int i = 0; i < m.rows(); ++i) { @@ -160,7 +186,7 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm s << fmt.rowSeparator; } s << fmt.matSuffix; - if(fmt.precision != -1) s.precision(old_precision); + if(explicit_precision) s.precision(old_precision); return s; } diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 16091caf0..1570f01e0 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -25,7 +25,7 @@ #ifndef EIGEN_MATHFUNCTIONS_H #define EIGEN_MATHFUNCTIONS_H -template typename NumTraits::Real epsilon() +template inline typename NumTraits::Real epsilon() { return std::numeric_limits::Real>::epsilon(); } From a5f820b873d61868c6833a004c56b8ba38c40d5d Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 20:03:14 -0400 Subject: [PATCH 29/58] forgot to update this --- Eigen/src/Core/util/Macros.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index f3f2b9066..dc3b3ee0a 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -244,7 +244,7 @@ using Eigen::ei_cos; // format used in Eigen's documentation // needed to define it here as escaping characters in CMake add_definition's argument seems very problematic. -#define EIGEN_DOCS_IO_FORMAT IOFormat(3, AlignCols, " ", "\n", "", "") +#define EIGEN_DOCS_IO_FORMAT IOFormat(3, 0, " ", "\n", "", "") #define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \ using Base::operator =; \ From fe4a86443f6f9b55c6da6e63c4227b0d9862f27f Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 20:16:04 -0400 Subject: [PATCH 30/58] fix warning --- Eigen/src/Core/IO.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h index 03e1af17d..a44802ed2 100644 --- a/Eigen/src/Core/IO.h +++ b/Eigen/src/Core/IO.h @@ -165,7 +165,7 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm width = std::max(width, int(sstr.str().length())); } } - std::streamsize old_precision; + std::streamsize old_precision = 0; if(explicit_precision) old_precision = s.precision(explicit_precision); s << fmt.matPrefix; for(int i = 0; i < m.rows(); ++i) From 8372bd12bde594af997977c1387bb4d2f8972ceb Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 20:23:13 -0400 Subject: [PATCH 31/58] update snippet --- doc/snippets/IOFormat.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/snippets/IOFormat.cpp b/doc/snippets/IOFormat.cpp index 9db70980c..735f5dd85 100644 --- a/doc/snippets/IOFormat.cpp +++ b/doc/snippets/IOFormat.cpp @@ -1,11 +1,11 @@ std::string sep = "\n----------------------------------------\n"; -Matrix3f m1; +Matrix3d m1; m1 << 1.111111, 2, 3.33333, 4, 5, 6, 7, 8.888888, 9; -IOFormat CommaInitFmt(4, Raw, ", ", ", ", "", "", " << ", ";"); -IOFormat CleanFmt(4, AlignCols, ", ", "\n", "[", "]"); -IOFormat OctaveFmt(4, AlignCols, ", ", ";\n", "", "", "[", "]"); -IOFormat HeavyFmt(4, AlignCols, ", ", ";\n", "[", "]", "[", "]"); +IOFormat CommaInitFmt(StreamPrecision, DontAlignCols, ", ", ", ", "", "", " << ", ";"); +IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]"); +IOFormat OctaveFmt(StreamPrecision, 0, ", ", ";\n", "", "", "[", "]"); +IOFormat HeavyFmt(FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]"); std::cout << m1 << sep; std::cout << m1.format(CommaInitFmt) << sep; From 62748a0963efa37a1e5f9b4fe2132403f16bc8fd Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 19:07:01 -0400 Subject: [PATCH 32/58] new script to generate and upload the docs for a given branch needs cleanup by a better shell scripter!! --- scripts/eigen_gen_docs | 63 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 scripts/eigen_gen_docs diff --git a/scripts/eigen_gen_docs b/scripts/eigen_gen_docs new file mode 100644 index 000000000..bc514b5a8 --- /dev/null +++ b/scripts/eigen_gen_docs @@ -0,0 +1,63 @@ +#!/bin/sh + +# todo: check that at most one argument was passed. + +if [ $# -eq 1 ] +then + branch=$1 + destination=dox-$1 +elif [ $# -eq 0 ] +then + branch=default + destination=dox-devel +fi + + +echo branch: $branch +echo destination: $destination + +# todo: push the `pwd` to restore it at the end + +rm -rf eigen_gen_docs_work_directory +mkdir eigen_gen_docs_work_directory + +# todo: is it really good practice to cd all the time? + +hg archive -r $branch eigen_gen_docs_work_directory/eigen2 + +cd eigen_gen_docs_work_directory + +#todo: check that the work directory was successfully created (exit if already existing) and entered + +mkdir build +cd build +cmake ../eigen2 +make -j3 doc #todo: n+1 where n = number of cpus + +#todo: check that make succeeded, is there a return code or something to check? + +#todo: check that there exists a doc subdir +cd doc + +#todo: check that there exists a html subdir +tar cfjv html.tar.bz2 html/ + +echo "put html.tar.bz2" > sftp_batchfile + +echo "uploading the dox archive" + +sftp ssh.tuxfamily.org < sftp_batchfile + +echo "cd eigen/eigen.tuxfamily.org-web/htdocs +pwd +mv ../../../html.tar.bz2 . +echo uncompressing the dox archive +tar xfjv html.tar.bz2 +echo removing old dox +rm -rf $destination +mv html $destination +echo giving write permissions to the group +chmod -R g+w $destination +" > ssh_batchfile + +ssh ssh.tuxfamily.org < ssh_batchfile \ No newline at end of file From 5fe0c308114eb24b2d5b7372d3c2c14c6f515b64 Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Sat, 15 Aug 2009 03:36:37 +0200 Subject: [PATCH 33/58] new script that update from mercurial, make the doc, and upload the result to tuxfamily.org --- scripts/eigen_gen_docs | 73 ++++++++++-------------------------------- 1 file changed, 17 insertions(+), 56 deletions(-) diff --git a/scripts/eigen_gen_docs b/scripts/eigen_gen_docs index bc514b5a8..3cdacc1a8 100644 --- a/scripts/eigen_gen_docs +++ b/scripts/eigen_gen_docs @@ -1,63 +1,24 @@ #!/bin/sh -# todo: check that at most one argument was passed. +# configuration +USER='orzel' -if [ $# -eq 1 ] +# step 1 : update +hg pull -u || (echo "update failed"; exit 1) + +# step 2 : build +# todo if 'build is not there, create one: +#mkdir build +(cd build && cmake .. && make -j3 doc) || (echo "make failed"; exit 1) +#todo: n+1 where n = number of cpus + +#step 3 : upload +BRANCH=`hg branch` +if [ $BRANCH == "default" ] then - branch=$1 - destination=dox-$1 -elif [ $# -eq 0 ] -then - branch=default - destination=dox-devel + BRANCH='devel' fi +# (the '/' at the end of path are very important, see rsync documentation) +rsync -az build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-$BRANCH/ || (echo "upload failed"; exit 1) -echo branch: $branch -echo destination: $destination - -# todo: push the `pwd` to restore it at the end - -rm -rf eigen_gen_docs_work_directory -mkdir eigen_gen_docs_work_directory - -# todo: is it really good practice to cd all the time? - -hg archive -r $branch eigen_gen_docs_work_directory/eigen2 - -cd eigen_gen_docs_work_directory - -#todo: check that the work directory was successfully created (exit if already existing) and entered - -mkdir build -cd build -cmake ../eigen2 -make -j3 doc #todo: n+1 where n = number of cpus - -#todo: check that make succeeded, is there a return code or something to check? - -#todo: check that there exists a doc subdir -cd doc - -#todo: check that there exists a html subdir -tar cfjv html.tar.bz2 html/ - -echo "put html.tar.bz2" > sftp_batchfile - -echo "uploading the dox archive" - -sftp ssh.tuxfamily.org < sftp_batchfile - -echo "cd eigen/eigen.tuxfamily.org-web/htdocs -pwd -mv ../../../html.tar.bz2 . -echo uncompressing the dox archive -tar xfjv html.tar.bz2 -echo removing old dox -rm -rf $destination -mv html $destination -echo giving write permissions to the group -chmod -R g+w $destination -" > ssh_batchfile - -ssh ssh.tuxfamily.org < ssh_batchfile \ No newline at end of file From 2f027b0d2ccfb396160fcbb8d21fd457d40f12e4 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 14 Aug 2009 21:58:41 -0400 Subject: [PATCH 34/58] only append the changeset to the version if we're in the default branch --- CMakeLists.txt | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d9a999c4..5d06f616a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,7 @@ project(Eigen) +cmake_minimum_required(VERSION 2.6.2) + # automatically parse the version number file(READ "${CMAKE_SOURCE_DIR}/Eigen/src/Core/util/Macros.h" _eigen2_version_header LIMIT 5000 OFFSET 1000) string(REGEX MATCH "define *EIGEN_WORLD_VERSION ([0-9]*)" _eigen2_world_version_match "${_eigen2_version_header}") @@ -10,21 +12,22 @@ string(REGEX MATCH "define *EIGEN_MINOR_VERSION ([0-9]*)" _eigen2_minor_version_ set(EIGEN2_MINOR_VERSION "${CMAKE_MATCH_1}") set(EIGEN_VERSION_NUMBER ${EIGEN2_WORLD_VERSION}.${EIGEN2_MAJOR_VERSION}.${EIGEN2_MINOR_VERSION}) -# if the mercurial program is absent, this will leave the EIGEN_HG_REVISION string empty, +# if the mercurial program is absent, this will leave the EIGEN_HG_CHANGESET string empty, # but won't stop CMake. execute_process(COMMAND hg tip -R ${CMAKE_SOURCE_DIR} OUTPUT_VARIABLE EIGEN_HGTIP_OUTPUT) +execute_process(COMMAND hg branch -R ${CMAKE_SOURCE_DIR} OUTPUT_VARIABLE EIGEN_BRANCH_OUTPUT) -# extract the mercurial revision number from the hg tip output -string(REGEX MATCH "^changeset: *[0-9]*:([0-9;a-f]+).*" EIGEN_HG_REVISION_MATCH "${EIGEN_HGTIP_OUTPUT}") -set(EIGEN_HG_REVISION "${CMAKE_MATCH_1}") - -if(EIGEN_HG_REVISION) - set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER} (mercurial revision ${EIGEN_HG_REVISION})") -else(EIGEN_HG_REVISION) +# if this is the default (aka development) branch, extract the mercurial changeset number from the hg tip output... +if(EIGEN_BRANCH_OUTPUT MATCHES "default") +string(REGEX MATCH "^changeset: *[0-9]*:([0-9;a-f]+).*" EIGEN_HG_CHANGESET_MATCH "${EIGEN_HGTIP_OUTPUT}") +set(EIGEN_HG_CHANGESET "${CMAKE_MATCH_1}") +endif(EIGEN_BRANCH_OUTPUT MATCHES "default") +#...and show it next to the version number +if(EIGEN_HG_CHANGESET) + set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER} (mercurial changeset ${EIGEN_HG_CHANGESET})") +else(EIGEN_HG_CHANGESET) set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}") -endif(EIGEN_HG_REVISION) - -cmake_minimum_required(VERSION 2.6.2) +endif(EIGEN_HG_CHANGESET) include(CheckCXXCompilerFlag) From 846e8b49ba15563729ef61b0645c7546d5f62db5 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 10:18:05 +0200 Subject: [PATCH 35/58] fix compilation of unit tests --- test/main.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/main.h b/test/main.h index 3947451cc..9ad014b64 100644 --- a/test/main.h +++ b/test/main.h @@ -45,7 +45,7 @@ namespace Eigen #define EI_PP_MAKE_STRING2(S) #S #define EI_PP_MAKE_STRING(S) EI_PP_MAKE_STRING2(S) -#define EIGEN_DEFAULT_IO_FORMAT IOFormat(4, AlignCols, " ", "\n", "", "", "", "") +#define EIGEN_DEFAULT_IO_FORMAT IOFormat(4, 0, " ", "\n", "", "", "", "") #ifndef EIGEN_NO_ASSERTION_CHECKING From d2becb9612c1cd2c8c45e1177b04d5543be8fb9d Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 10:19:16 +0200 Subject: [PATCH 36/58] add a "rot" benchmark in BTL --- bench/btl/actions/action_rot.hh | 116 ++++++++++++++++++++++ bench/btl/actions/basic_actions.hh | 1 + bench/btl/data/action_settings.txt | 3 +- bench/btl/data/go_mean | 1 + bench/btl/libs/C_BLAS/C_BLAS_interface.hh | 8 ++ bench/btl/libs/C_BLAS/main.cpp | 1 + bench/btl/libs/eigen2/eigen2_interface.hh | 4 + bench/btl/libs/eigen2/main_linear.cpp | 1 + 8 files changed, 134 insertions(+), 1 deletion(-) create mode 100644 bench/btl/actions/action_rot.hh diff --git a/bench/btl/actions/action_rot.hh b/bench/btl/actions/action_rot.hh new file mode 100644 index 000000000..df822a6d6 --- /dev/null +++ b/bench/btl/actions/action_rot.hh @@ -0,0 +1,116 @@ + +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License +// as published by the Free Software Foundation; either version 2 +// of the License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +// +#ifndef ACTION_ROT +#define ACTION_ROT +#include "utilities.h" +#include "STL_interface.hh" +#include +#include "init/init_function.hh" +#include "init/init_vector.hh" +#include "init/init_matrix.hh" + +using namespace std; + +template +class Action_rot { + +public : + + // Ctor + BTL_DONT_INLINE Action_rot( int size ):_size(size) + { + MESSAGE("Action_rot Ctor"); + + // STL matrix and vector initialization + typename Interface::stl_matrix tmp; + init_vector(A_stl,_size); + init_vector(B_stl,_size); + + // generic matrix and vector initialization + Interface::vector_from_stl(A_ref,A_stl); + Interface::vector_from_stl(A,A_stl); + Interface::vector_from_stl(B_ref,B_stl); + Interface::vector_from_stl(B,B_stl); + } + + // invalidate copy ctor + Action_rot( const Action_rot & ) + { + INFOS("illegal call to Action_rot Copy Ctor"); + exit(1); + } + + // Dtor + BTL_DONT_INLINE ~Action_rot( void ){ + MESSAGE("Action_rot Dtor"); + Interface::free_vector(A); + Interface::free_vector(B); + Interface::free_vector(A_ref); + Interface::free_vector(B_ref); + } + + // action name + static inline std::string name( void ) + { + return "rot_" + Interface::name(); + } + + double nb_op_base( void ){ + return 6.0*_size; + } + + BTL_DONT_INLINE void initialize( void ){ + Interface::copy_vector(A_ref,A,_size); + Interface::copy_vector(B_ref,B,_size); + } + + BTL_DONT_INLINE void calculate( void ) { + BTL_ASM_COMMENT("#begin rot"); + Interface::rot(A,B,0.5,0.6,_size); + BTL_ASM_COMMENT("end rot"); + } + + BTL_DONT_INLINE void check_result( void ){ + // calculation check +// Interface::vector_to_stl(X,resu_stl); + +// STL_interface::rot(A_stl,B_stl,X_stl,_size); + +// typename Interface::real_type error= +// STL_interface::norm_diff(X_stl,resu_stl); + +// if (error>1.e-3){ +// INFOS("WRONG CALCULATION...residual=" << error); +// exit(0); +// } + + } + +private : + + typename Interface::stl_vector A_stl; + typename Interface::stl_vector B_stl; + + typename Interface::gene_vector A_ref; + typename Interface::gene_vector B_ref; + + typename Interface::gene_vector A; + typename Interface::gene_vector B; + + int _size; +}; + + +#endif diff --git a/bench/btl/actions/basic_actions.hh b/bench/btl/actions/basic_actions.hh index 06a382be6..beaeea8c3 100644 --- a/bench/btl/actions/basic_actions.hh +++ b/bench/btl/actions/basic_actions.hh @@ -15,6 +15,7 @@ // #include "action_symm.hh" #include "action_syr2.hh" #include "action_ger.hh" +#include "action_rot.hh" // #include "action_lu_solve.hh" diff --git a/bench/btl/data/action_settings.txt b/bench/btl/data/action_settings.txt index e2d931ddc..9bee1651c 100644 --- a/bench/btl/data/action_settings.txt +++ b/bench/btl/data/action_settings.txt @@ -14,4 +14,5 @@ tridiagonalization ; "{/*1.5 Tridiagonalization}" ; "matrix size" ; 4:3000 hessenberg ; "{/*1.5 Hessenberg decomposition}" ; "matrix size" ; 4:3000 symv ; "{/*1.5 symmetric matrix vector product}" ; "matrix size" ; 4:3000 syr2 ; "{/*1.5 symmretric rank-2 update (A += u^T v + u v^T)}" ; "matrix size" ; 4:3000 -ger ; "{/*1.5 general rank-1 update (A += u v^T)}" ; "matrix size" ; 4:3000 \ No newline at end of file +ger ; "{/*1.5 general rank-1 update (A += u v^T)}" ; "matrix size" ; 4:3000 +rot ; "{/*1.5 apply rotation in the plane}" ; "vector size" ; 4:1000000 \ No newline at end of file diff --git a/bench/btl/data/go_mean b/bench/btl/data/go_mean index 5ff3a3b98..4e4fd295a 100755 --- a/bench/btl/data/go_mean +++ b/bench/btl/data/go_mean @@ -48,6 +48,7 @@ source mk_mean_script.sh hessenberg $1 11 100 300 1000 $mode $prefix source mk_mean_script.sh symv $1 11 50 300 1000 $mode $prefix source mk_mean_script.sh syr2 $1 11 50 300 1000 $mode $prefix source mk_mean_script.sh ger $1 11 50 300 1000 $mode $prefix +source mk_mean_script.sh rot $1 11 2500 100000 250000 $mode $prefix fi diff --git a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh index 725f35944..f07df1033 100644 --- a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh +++ b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh @@ -179,6 +179,14 @@ public : #endif } + static inline void rot(gene_vector & A, gene_vector & B, float c, float s, int N){ + #ifdef PUREBLAS + srot_(&N,A,&intone,B,&intone,&c,&s); + #else + cblas_srot(N,A,1,B,1,c,s); + #endif + } + static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ #ifdef PUREBLAS sgemv_(&trans,&N,&N,&fone,A,&N,B,&intone,&fzero,X,&intone); diff --git a/bench/btl/libs/C_BLAS/main.cpp b/bench/btl/libs/C_BLAS/main.cpp index 6fd75812c..fdffdc979 100644 --- a/bench/btl/libs/C_BLAS/main.cpp +++ b/bench/btl/libs/C_BLAS/main.cpp @@ -45,6 +45,7 @@ int main() bench > >(MIN_MV,MAX_MV,NB_POINT); bench > >(MIN_MV,MAX_MV,NB_POINT); + bench > >(MIN_AXPY,MAX_AXPY,NB_POINT); bench > >(MIN_MM,MAX_MM,NB_POINT); bench > >(MIN_MM,MAX_MM,NB_POINT); diff --git a/bench/btl/libs/eigen2/eigen2_interface.hh b/bench/btl/libs/eigen2/eigen2_interface.hh index 6a1bc5d61..1a5f89834 100644 --- a/bench/btl/libs/eigen2/eigen2_interface.hh +++ b/bench/btl/libs/eigen2/eigen2_interface.hh @@ -168,6 +168,10 @@ public : A.col(j) += X * Y[j]; } + static EIGEN_DONT_INLINE void rot(gene_vector & A, gene_vector & B, real c, real s, int N){ + ei_apply_rotation_in_the_plane(A, B, c, s); + } + static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ X = (A.transpose()*B).lazy(); } diff --git a/bench/btl/libs/eigen2/main_linear.cpp b/bench/btl/libs/eigen2/main_linear.cpp index e79927b0f..e73fe2043 100644 --- a/bench/btl/libs/eigen2/main_linear.cpp +++ b/bench/btl/libs/eigen2/main_linear.cpp @@ -27,6 +27,7 @@ int main() bench > >(MIN_AXPY,MAX_AXPY,NB_POINT); bench > >(MIN_AXPY,MAX_AXPY,NB_POINT); + bench > >(MIN_AXPY,MAX_AXPY,NB_POINT); return 0; } From 109a4f650bb69f1565f9562b384d711e36c56b26 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 10:20:01 +0200 Subject: [PATCH 37/58] fix a couple of warnings --- Eigen/src/Core/MapBase.h | 4 ++-- Eigen/src/Core/StableNorm.h | 12 ++++++------ Eigen/src/Core/TriangularMatrix.h | 4 ++-- test/product_trsm.cpp | 6 +++--- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index 0151cc67f..51e041db6 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -66,11 +66,11 @@ template class MapBase inline const Scalar* data() const { return m_data; } template struct force_aligned_impl { - AlignedDerivedType static run(MapBase& a) { return a.derived(); } + static AlignedDerivedType run(MapBase& a) { return a.derived(); } }; template struct force_aligned_impl { - AlignedDerivedType static run(MapBase& a) { return a.derived()._convertToForceAligned(); } + static AlignedDerivedType run(MapBase& a) { return a.derived()._convertToForceAligned(); } }; /** \returns an expression equivalent to \c *this but having the \c PacketAccess constant diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 22809633d..6056cafab 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -115,20 +115,20 @@ MatrixBase::blueNorm() const ei_assert(false && "the algorithm cannot be guaranteed on this computer"); } iexp = -((1-iemin)/2); - b1 = std::pow(double(ibeta),iexp); // lower boundary of midrange + b1 = Scalar(std::pow(double(ibeta),iexp)); // lower boundary of midrange iexp = (iemax + 1 - it)/2; - b2 = std::pow(double(ibeta),iexp); // upper boundary of midrange + b2 = Scalar(std::pow(double(ibeta),iexp)); // upper boundary of midrange iexp = (2-iemin)/2; - s1m = std::pow(double(ibeta),iexp); // scaling factor for lower range + s1m = Scalar(std::pow(double(ibeta),iexp)); // scaling factor for lower range iexp = - ((iemax+it)/2); - s2m = std::pow(double(ibeta),iexp); // scaling factor for upper range + s2m = Scalar(std::pow(double(ibeta),iexp)); // scaling factor for upper range overfl = rbig*s2m; // overfow boundary for abig - eps = std::pow(double(ibeta), 1-it); + eps = Scalar(std::pow(double(ibeta), 1-it)); relerr = ei_sqrt(eps); // tolerance for neglecting asml abig = 1.0/eps - 1.0; - if (RealScalar(nbig)>abig) nmax = abig; // largest safe n + if (RealScalar(nbig)>abig) nmax = int(abig); // largest safe n else nmax = nbig; } int n = size(); diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index c262ea7a7..b0362f20c 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -234,14 +234,14 @@ template class TriangularView inline TriangularView,TransposeMode> adjoint() { return m_matrix.adjoint().nestByValue(); } /** \sa MatrixBase::adjoint() const */ - const inline TriangularView,TransposeMode> adjoint() const + inline const TriangularView,TransposeMode> adjoint() const { return m_matrix.adjoint().nestByValue(); } /** \sa MatrixBase::transpose() */ inline TriangularView >,TransposeMode> transpose() { return m_matrix.transpose().nestByValue(); } /** \sa MatrixBase::transpose() const */ - const inline TriangularView >,TransposeMode> transpose() const + inline const TriangularView >,TransposeMode> transpose() const { return m_matrix.transpose().nestByValue(); } PlainMatrixType toDense() const diff --git a/test/product_trsm.cpp b/test/product_trsm.cpp index 4f0fd15be..9f372df91 100644 --- a/test/product_trsm.cpp +++ b/test/product_trsm.cpp @@ -25,9 +25,9 @@ #include "main.h" #define VERIFY_TRSM(TRI,XB) { \ - XB.setRandom(); ref = XB; \ - TRI.template solveInPlace(XB); \ - VERIFY_IS_APPROX(TRI.toDense() * XB, ref); \ + (XB).setRandom(); ref = (XB); \ + (TRI).solveInPlace(XB); \ + VERIFY_IS_APPROX((TRI).toDense() * (XB), ref); \ } template void trsm(int size,int cols) From bff4238d154c1bbb37abd5f5df1bf04f0117e2cc Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 10:24:27 +0200 Subject: [PATCH 38/58] fix setFromTwoVectors because of the change in sorting of the the singular values --- Eigen/src/Geometry/Quaternion.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 2f9f97807..3d59533f4 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -373,7 +373,7 @@ inline Quaternion& Quaternion::setFromTwoVectors(const MatrixBas c = std::max(c,-1); Matrix m; m << v0.transpose(), v1.transpose(); SVD > svd(m); - Vector3 axis = svd.matrixV().col(2); + Vector3 axis = svd.matrixV().col(0); Scalar w2 = (Scalar(1)+c)*Scalar(0.5); this->w() = ei_sqrt(w2); From 0da31a6e1dfa06ec87b967a4124d03c3cf3a4389 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 10:55:11 +0200 Subject: [PATCH 39/58] bugfix and compilation fix in ProductBase --- Eigen/src/Core/ProductBase.h | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index 0da046b1e..2090a0e49 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -115,7 +115,7 @@ class ProductBase : public MatrixBase { PlainMatrixType res(rows(), cols()); res.setZero(); - evalTo(res); + derived().evalTo(res); return res; } @@ -148,12 +148,14 @@ class ScaledProduct; // functions of ProductBase, because, otherwise we would have to // define all overloads defined in MatrixBase. Furthermore, Using // "using Base::operator*" would not work with MSVC. -template -const ScaledProduct operator*(const ProductBase& prod, typename Derived::Scalar x) +// +// Also note that here we accept any type which can be converted to Derived::Scalar. +template +const ScaledProduct operator*(const ProductBase& prod, Scalar x) { return ScaledProduct(prod.derived(), x); } -template -const ScaledProduct operator*(typename Derived::Scalar x,const ProductBase& prod) +template +const ScaledProduct operator*(Scalar x,const ProductBase& prod) { return ScaledProduct(prod.derived(), x); } template @@ -176,7 +178,7 @@ class ScaledProduct typedef typename Base::Scalar Scalar; // EIGEN_PRODUCT_PUBLIC_INTERFACE(ScaledProduct) - ScaledProduct(const NestedProduct& prod, Scalar& x) + ScaledProduct(const NestedProduct& prod, Scalar x) : Base(prod.lhs(),prod.rhs()), m_prod(prod), m_alpha(x) {} template From 7b60713e87ba24d3b0fa8f486bba8665151e9f31 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 11:52:50 +0200 Subject: [PATCH 40/58] my previous fix was not very good --- Eigen/src/Core/ProductBase.h | 26 +++++++++++++++++++++----- Eigen/src/Core/util/Meta.h | 7 +------ 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index 2090a0e49..b5ed5ae8f 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -149,15 +149,31 @@ class ScaledProduct; // define all overloads defined in MatrixBase. Furthermore, Using // "using Base::operator*" would not work with MSVC. // -// Also note that here we accept any type which can be converted to Derived::Scalar. -template -const ScaledProduct operator*(const ProductBase& prod, Scalar x) +// Also note that here we accept any compatible scalar types +template +const ScaledProduct +operator*(const ProductBase& prod, typename Derived::Scalar x) { return ScaledProduct(prod.derived(), x); } -template -const ScaledProduct operator*(Scalar x,const ProductBase& prod) +template +typename ei_enable_if::ret, + const ScaledProduct >::type +operator*(const ProductBase& prod, typename Derived::RealScalar x) { return ScaledProduct(prod.derived(), x); } + +template +const ScaledProduct +operator*(typename Derived::Scalar x,const ProductBase& prod) +{ return ScaledProduct(prod.derived(), x); } + +template +typename ei_enable_if::ret, + const ScaledProduct >::type +operator*(typename Derived::RealScalar x,const ProductBase& prod) +{ return ScaledProduct(prod.derived(), x); } + + template struct ei_traits > : ei_traits, diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 029d966fd..3a960bea6 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -155,12 +155,7 @@ template class ei_meta_sqrt { public: enum { ret = (SupX*SupX <= Y) ? SupX : InfX }; }; /** \internal determines whether the product of two numeric types is allowed and what the return type is */ -template struct ei_scalar_product_traits -{ - // dummy general case where T and U aren't compatible -- not allowed anyway but we catch it elsewhere - //enum { Cost = NumTraits::MulCost }; - typedef T ReturnType; -}; +template struct ei_scalar_product_traits; template struct ei_scalar_product_traits { From 13a8956188d1c8441c01a6218dc24c745ff0b210 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 13:12:50 +0200 Subject: [PATCH 41/58] bugfix in inner-product specialization, compilation fix in stable norm, optimize apply householder --- Eigen/src/Core/Product.h | 2 +- Eigen/src/Core/StableNorm.h | 10 +++++----- Eigen/src/Householder/Householder.h | 6 ++++-- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 610d5c84a..fba4241b1 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -156,7 +156,7 @@ class GeneralProduct template void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(dst.rows()==1 && dst.cols()==1); - dst.coeffRef(0,0) += alpha * (m_lhs.cwise()*m_rhs).sum(); + dst.coeffRef(0,0) += alpha * (m_lhs.transpose().cwise()*m_rhs).sum(); } }; diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 6056cafab..77fe79782 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -115,17 +115,17 @@ MatrixBase::blueNorm() const ei_assert(false && "the algorithm cannot be guaranteed on this computer"); } iexp = -((1-iemin)/2); - b1 = Scalar(std::pow(double(ibeta),iexp)); // lower boundary of midrange + b1 = RealScalar(std::pow(double(ibeta),iexp)); // lower boundary of midrange iexp = (iemax + 1 - it)/2; - b2 = Scalar(std::pow(double(ibeta),iexp)); // upper boundary of midrange + b2 = RealScalar(std::pow(double(ibeta),iexp)); // upper boundary of midrange iexp = (2-iemin)/2; - s1m = Scalar(std::pow(double(ibeta),iexp)); // scaling factor for lower range + s1m = RealScalar(std::pow(double(ibeta),iexp)); // scaling factor for lower range iexp = - ((iemax+it)/2); - s2m = Scalar(std::pow(double(ibeta),iexp)); // scaling factor for upper range + s2m = RealScalar(std::pow(double(ibeta),iexp)); // scaling factor for upper range overfl = rbig*s2m; // overfow boundary for abig - eps = Scalar(std::pow(double(ibeta), 1-it)); + eps = RealScalar(std::pow(double(ibeta), 1-it)); relerr = ei_sqrt(eps); // tolerance for neglecting asml abig = 1.0/eps - 1.0; if (RealScalar(nbig)>abig) nmax = int(abig); // largest safe n diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 2a4b1395c..42635cb18 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -73,7 +73,8 @@ void MatrixBase::applyHouseholderOnTheLeft( { Matrix tmp(cols()); Block bottom(derived(), 1, 0, rows()-1, cols()); - tmp = row(0) + essential.adjoint() * bottom; + tmp = (essential.adjoint() * bottom).lazy(); + tmp += row(0); row(0) -= beta * tmp; bottom -= beta * essential * tmp; } @@ -86,7 +87,8 @@ void MatrixBase::applyHouseholderOnTheRight( { Matrix tmp(rows()); Block right(derived(), 0, 1, rows(), cols()-1); - tmp = col(0) + right * essential.conjugate(); + tmp = (right * essential.conjugate()).lazy(); + tmp += col(0); col(0) -= beta * tmp; right -= beta * tmp * essential.transpose(); } From 50c703f0c750ebbf1aaec785dcddeb0fb763877f Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 18:35:51 +0200 Subject: [PATCH 42/58] As proposed on the list: - rename EvalBeforeAssignBit to MayAliasBit - make .lazy() remove the MayAliasBit only, and mark it as deprecated - add a NoAlias pseudo expression, and MatrixBase::noalias() function Todo: - we have to decide whether += and -= assume no aliasing by default ? - once we agree on the API: update the Sparse module and the unit tests respectively. --- Eigen/Core | 1 + Eigen/src/Core/Assign.h | 2 +- Eigen/src/Core/Flagged.h | 15 ++-- Eigen/src/Core/MapBase.h | 4 +- Eigen/src/Core/Matrix.h | 2 +- Eigen/src/Core/MatrixBase.h | 13 ++-- Eigen/src/Core/NoAlias.h | 94 +++++++++++++++++++++++ Eigen/src/Core/ProductBase.h | 8 +- Eigen/src/Core/TriangularMatrix.h | 6 +- Eigen/src/Core/products/GeneralUnrolled.h | 2 +- Eigen/src/Core/util/Constants.h | 8 +- Eigen/src/Core/util/ForwardDeclarations.h | 1 + test/product_notemporary.cpp | 8 +- test/submatrices.cpp | 3 +- 14 files changed, 138 insertions(+), 29 deletions(-) create mode 100644 Eigen/src/Core/NoAlias.h diff --git a/Eigen/Core b/Eigen/Core index 28ed09035..ff56c6ca2 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -154,6 +154,7 @@ namespace Eigen { #include "src/Core/NestByValue.h" #include "src/Core/ReturnByValue.h" #include "src/Core/Flagged.h" +#include "src/Core/NoAlias.h" #include "src/Core/Matrix.h" #include "src/Core/Cwise.h" #include "src/Core/CwiseBinaryOp.h" diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 4bd1046a7..b0e69224e 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -434,7 +434,7 @@ EIGEN_STRONG_INLINE Derived& MatrixBase } template clas ExpressionTypeNested m_matrix; }; -/** \returns an expression of *this with added flags +/** \deprecated it is only used by lazy() which is deprecated * - * \addexample MarkExample \label How to mark a triangular matrix as triangular + * \returns an expression of *this with added flags * * Example: \include MatrixBase_marked.cpp * Output: \verbinclude MatrixBase_marked.out @@ -128,8 +130,9 @@ MatrixBase::marked() const return derived(); } -/** \returns an expression of *this with the following flags removed: - * EvalBeforeNestingBit and EvalBeforeAssigningBit. +/** \deprecated use MatrixBase::noalias() + * + * \returns an expression of *this with the MayAliasBit flag removed. * * Example: \include MatrixBase_lazy.cpp * Output: \verbinclude MatrixBase_lazy.out @@ -137,7 +140,7 @@ MatrixBase::marked() const * \sa class Flagged, marked() */ template -inline const Flagged +inline const Flagged MatrixBase::lazy() const { return derived(); diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index 51e041db6..b9e9374be 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -179,11 +179,11 @@ template class MapBase // explicitly add these two overloads. // Maybe there exists a better solution though. template - Derived& operator+=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) + Derived& operator+=(const Flagged, 0, MayAliasBit>& other) { return Base::operator+=(other); } template - Derived& operator-=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) + Derived& operator-=(const Flagged, 0, MayAliasBit>& other) { return Base::operator-=(other); } template diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index f58424ba2..d129535b6 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -577,7 +577,7 @@ class Matrix template EIGEN_STRONG_INLINE Matrix& _set(const MatrixBase& other) { - _set_selector(other.derived(), typename ei_meta_if(int(OtherDerived::Flags) & EvalBeforeAssigningBit), ei_meta_true, ei_meta_false>::ret()); + _set_selector(other.derived(), typename ei_meta_if(int(OtherDerived::Flags) & MayAliasBit), ei_meta_true, ei_meta_false>::ret()); return *this; } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 52d5f680c..a554de8e3 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -326,9 +326,10 @@ template class MatrixBase template Derived& lazyAssign(const MatrixBase& other); - /** Overloaded for cache friendly product evaluation */ + /** \deprecated because .lazy() is deprecated + * Overloaded for cache friendly product evaluation */ template - Derived& lazyAssign(const Flagged& other) + Derived& lazyAssign(const Flagged& other) { return lazyAssign(other._expression()); } template @@ -336,11 +337,11 @@ template class MatrixBase template Derived& operator+=(const Flagged, 0, - EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); + MayAliasBit>& other); template Derived& operator-=(const Flagged, 0, - EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); + MayAliasBit>& other); #endif // not EIGEN_PARSED_BY_DOXYGEN CommaInitializer operator<< (const Scalar& s); @@ -614,7 +615,9 @@ template class MatrixBase template const Flagged marked() const; - const Flagged lazy() const; + const Flagged lazy() const; + + NoAlias noalias(); /** \returns number of elements to skip to pass from one row (resp. column) to another * for a row-major (resp. column-major) matrix. diff --git a/Eigen/src/Core/NoAlias.h b/Eigen/src/Core/NoAlias.h new file mode 100644 index 000000000..5c2b6e3ca --- /dev/null +++ b/Eigen/src/Core/NoAlias.h @@ -0,0 +1,94 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_NOALIAS_H +#define EIGEN_NOALIAS_H + +/** \class NoAlias + * + * \brief Pseudo expression providing an operator = assuming no aliasing + * + * \param ExpressionType the type of the object on which to do the lazy assignment + * + * This class represents an expression with a special assignment operator (operator=) + * assuming no aliasing between the target expression and the source expression. + * It is the return type of MatrixBase::noalias() + * and most of the time this is the only way it is used. + * + * \sa MatrixBase::noalias() + */ +template +class NoAlias +{ + public: + NoAlias(ExpressionType& expression) : m_expression(expression) {} + + /** Behaves like MatrixBase::lazyAssign() */ + template + ExpressionType& operator=(const MatrixBase& other) + { + return m_expression.lazyAssign(other.derived()); + } + + // TODO could be removed if we decide that += is noalias by default + template + ExpressionType& operator+=(const MatrixBase& other) + { + return m_expression.lazyAssign(m_expression + other.derived()); + } + + // TODO could be removed if we decide that += is noalias by default + template + ExpressionType& operator-=(const MatrixBase& other) + { + return m_expression.lazyAssign(m_expression - other.derived()); + } + + // TODO could be removed if we decide that += is noalias by default + template + ExpressionType& operator+=(const ProductBase& other) + { other.derived().addTo(m_expression); return m_expression; } + + // TODO could be removed if we decide that += is noalias by default + template + ExpressionType& operator-=(const ProductBase& other) + { other.derived().subTo(m_expression); return m_expression; } + + protected: + ExpressionType& m_expression; +}; + + +/** \returns a pseudo expression of \c *this with an operator= assuming + * no aliasing between \c *this and the source expression + * + * \sa class NoAlias + */ +template +NoAlias MatrixBase::noalias() +{ + return derived(); +} + +#endif // EIGEN_NOALIAS_H diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index b5ed5ae8f..867aa7a15 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -39,7 +39,7 @@ struct ei_traits > ColsAtCompileTime = ei_traits::ColsAtCompileTime, MaxRowsAtCompileTime = ei_traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = ei_traits::MaxColsAtCompileTime, - Flags = EvalBeforeNestingBit | EvalBeforeAssigningBit, + Flags = EvalBeforeNestingBit | MayAliasBit, CoeffReadCost = 0 // FIXME why is it needed ? }; }; @@ -119,7 +119,7 @@ class ProductBase : public MatrixBase return res; } - const Flagged lazy() const + const Flagged lazy() const { return *this; } @@ -228,7 +228,7 @@ Derived& MatrixBase::lazyAssign(const ProductBase template Derived& MatrixBase::operator+=(const Flagged, 0, - EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) + MayAliasBit>& other) { other._expression().derived().addTo(derived()); return derived(); } @@ -238,7 +238,7 @@ Derived& MatrixBase::operator+=(const Flagged template Derived& MatrixBase::operator-=(const Flagged, 0, - EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) + MayAliasBit>& other) { other._expression().derived().subTo(derived()); return derived(); } diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index b0362f20c..cfaf89733 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -477,7 +477,7 @@ template inline TriangularView& TriangularView::operator=(const MatrixBase& other) { - if(OtherDerived::Flags & EvalBeforeAssigningBit) + if(OtherDerived::Flags & MayAliasBit) { typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); other_evaluated.template triangularView().lazyAssign(other.derived()); @@ -512,7 +512,7 @@ inline TriangularView& TriangularView::operator=(const TriangularBase& other) { ei_assert(Mode == OtherDerived::Mode); - if(ei_traits::Flags & EvalBeforeAssigningBit) + if(ei_traits::Flags & MayAliasBit) { typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); other_evaluated.template triangularView().lazyAssign(other.derived()); @@ -548,7 +548,7 @@ template template void TriangularBase::evalToDense(MatrixBase &other) const { - if(ei_traits::Flags & EvalBeforeAssigningBit) + if(ei_traits::Flags & MayAliasBit) { typename Derived::PlainMatrixType other_evaluated(rows(), cols()); evalToDenseLazy(other_evaluated); diff --git a/Eigen/src/Core/products/GeneralUnrolled.h b/Eigen/src/Core/products/GeneralUnrolled.h index 7241976a8..2961d9452 100644 --- a/Eigen/src/Core/products/GeneralUnrolled.h +++ b/Eigen/src/Core/products/GeneralUnrolled.h @@ -76,7 +76,7 @@ struct ei_traits > RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) - | EvalBeforeAssigningBit + | MayAliasBit | EvalBeforeNestingBit | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0) | (LhsFlags & RhsFlags & AlignedBit), diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 21e7f62c3..1d55217f5 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -76,8 +76,10 @@ const unsigned int EvalBeforeNestingBit = 0x2; /** \ingroup flags * - * means the expression should be evaluated before any assignement */ -const unsigned int EvalBeforeAssigningBit = 0x4; + * Means the expression cannot be evaluated safely if the result alias one + * of the operands of the expression. Therefore it should be evaluated + * before any assignement. */ +const unsigned int MayAliasBit = 0x4; /** \ingroup flags * @@ -182,7 +184,7 @@ const unsigned int SparseBit = 0x1000; // list of flags that are inherited by default const unsigned int HereditaryBits = RowMajorBit | EvalBeforeNestingBit - | EvalBeforeAssigningBit + | MayAliasBit | SparseBit; // Possible values for the Mode parameter of part() diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index c9ce174d2..6f3e6a788 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -35,6 +35,7 @@ template class Matrix; template class Flagged; +template class NoAlias; template class NestByValue; template class SwapWrapper; template class Minor; diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index d5d996e49..f6e8f4101 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -71,15 +71,19 @@ template void product_notemporary(const MatrixType& m) VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()), 1); VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m1 * m2.adjoint(), 0); + // NOTE in this case the slow product is used: // FIXME: -// VERIFY_EVALUATION_COUNT( m3 = s1 * (m1 * m2.transpose()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * (m1 * m2.transpose()), 0); VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * m2.adjoint()).lazy(), 0); VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * (m1*s3+m2*s2).adjoint()).lazy(), 1); VERIFY_EVALUATION_COUNT( m3 = ((s1 * m1).adjoint() * s2 * m2).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 -= (s1 * (-m1*s3).adjoint() * (s2 * m2 * s3)).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3 += (s1 * (-m1*s3).adjoint() * (s2 * m2 * s3)).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() += s1 * (-m1*s3).adjoint() * (s2 * m2 * s3), 0); VERIFY_EVALUATION_COUNT( m3 -= (s1 * (m1.transpose() * m2)).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= s1 * (m1.transpose() * m2), 0); VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) += (-m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint()).lazy() ), 0); VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) -= (s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1)).lazy() ), 0); diff --git a/test/submatrices.cpp b/test/submatrices.cpp index fc7fd09ea..a819cadc2 100644 --- a/test/submatrices.cpp +++ b/test/submatrices.cpp @@ -86,7 +86,8 @@ template void submatrices(const MatrixType& m) //check row() and col() VERIFY_IS_APPROX(m1.col(c1).transpose(), m1.transpose().row(c1)); - VERIFY_IS_APPROX(m1.col(c1).dot(square.row(r1)), (square.lazy() * m1.conjugate()).eval()(r1,c1)); + // FIXME perhaps we should re-enable that without the .eval() + VERIFY_IS_APPROX(m1.col(c1).dot(square.row(r1)), (square * m1.conjugate()).eval()(r1,c1)); //check operator(), both constant and non-constant, on row() and col() m1.row(r1) += s1 * m1.row(r2); m1.col(c1) += s1 * m1.col(c2); From a3e6047c25a4cbc2153974e04fe124c5776a23c0 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sat, 15 Aug 2009 15:29:44 -0400 Subject: [PATCH 43/58] fix and improve docs --- Eigen/src/Array/Random.h | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Array/Random.h b/Eigen/src/Array/Random.h index 61c76bbd3..15cc6ae7c 100644 --- a/Eigen/src/Array/Random.h +++ b/Eigen/src/Array/Random.h @@ -35,13 +35,13 @@ struct ei_functor_traits > /** \array_module * - * \returns a random matrix (not an expression, the matrix is immediately evaluated). + * \returns a random matrix expression * * The parameters \a rows and \a cols are the number of rows and of columns of * the returned matrix. Must be compatible with this MatrixBase type. * * This variant is meant to be used for dynamic-size matrix types. For fixed-size types, - * it is redundant to pass \a rows and \a cols as arguments, so ei_random() should be used + * it is redundant to pass \a rows and \a cols as arguments, so Random() should be used * instead. * * \addexample RandomExample \label How to create a matrix with random coefficients @@ -49,6 +49,10 @@ struct ei_functor_traits > * Example: \include MatrixBase_random_int_int.cpp * Output: \verbinclude MatrixBase_random_int_int.out * + * This expression has the "evaluate before nesting" flag so that it will be evaluated into + * a temporary matrix whenever it is nested in a larger expression. This prevents unexpected + * behavior with expressions involving random matrices. + * * \sa MatrixBase::setRandom(), MatrixBase::Random(int), MatrixBase::Random() */ template @@ -60,7 +64,7 @@ MatrixBase::Random(int rows, int cols) /** \array_module * - * \returns a random vector (not an expression, the vector is immediately evaluated). + * \returns a random vector expression * * The parameter \a size is the size of the returned vector. * Must be compatible with this MatrixBase type. @@ -68,12 +72,16 @@ MatrixBase::Random(int rows, int cols) * \only_for_vectors * * This variant is meant to be used for dynamic-size vector types. For fixed-size types, - * it is redundant to pass \a size as argument, so ei_random() should be used + * it is redundant to pass \a size as argument, so Random() should be used * instead. * * Example: \include MatrixBase_random_int.cpp * Output: \verbinclude MatrixBase_random_int.out * + * This expression has the "evaluate before nesting" flag so that it will be evaluated into + * a temporary vector whenever it is nested in a larger expression. This prevents unexpected + * behavior with expressions involving random matrices. + * * \sa MatrixBase::setRandom(), MatrixBase::Random(int,int), MatrixBase::Random() */ template @@ -85,8 +93,7 @@ MatrixBase::Random(int size) /** \array_module * - * \returns a fixed-size random matrix or vector - * (not an expression, the matrix is immediately evaluated). + * \returns a fixed-size random matrix or vector expression * * This variant is only for fixed-size MatrixBase types. For dynamic-size types, you * need to use the variants taking size arguments. @@ -94,6 +101,10 @@ MatrixBase::Random(int size) * Example: \include MatrixBase_random.cpp * Output: \verbinclude MatrixBase_random.out * + * This expression has the "evaluate before nesting" flag so that it will be evaluated into + * a temporary matrix whenever it is nested in a larger expression. This prevents unexpected + * behavior with expressions involving random matrices. + * * \sa MatrixBase::setRandom(), MatrixBase::Random(int,int), MatrixBase::Random(int) */ template From 239ada95b7680c75f793086eaa35fe7ec1047204 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 22:19:29 +0200 Subject: [PATCH 44/58] add overloads of lazyAssign to detect common aliasing issue with transpose and adjoint --- Eigen/Core | 2 +- Eigen/src/Core/CwiseBinaryOp.h | 6 +++ Eigen/src/Core/MatrixBase.h | 15 ++++++ Eigen/src/Core/Transpose.h | 88 ++++++++++++++++++++++++++++++++++ test/adjoint.cpp | 14 +++++- 5 files changed, 122 insertions(+), 3 deletions(-) diff --git a/Eigen/Core b/Eigen/Core index ff56c6ca2..854f737d6 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -150,6 +150,7 @@ namespace Eigen { #include "src/Core/Assign.h" #endif +#include "src/Core/util/BlasUtil.h" #include "src/Core/MatrixStorage.h" #include "src/Core/NestByValue.h" #include "src/Core/ReturnByValue.h" @@ -178,7 +179,6 @@ namespace Eigen { #include "src/Core/IO.h" #include "src/Core/Swap.h" #include "src/Core/CommaInitializer.h" -#include "src/Core/util/BlasUtil.h" #include "src/Core/ProductBase.h" #include "src/Core/Product.h" #include "src/Core/TriangularMatrix.h" diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index bcd4ae25f..78de9fbce 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -85,6 +85,8 @@ class CwiseBinaryOp : ei_no_assignment_operator, EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) typedef typename ei_traits::LhsNested LhsNested; typedef typename ei_traits::RhsNested RhsNested; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::_RhsNested _RhsNested; EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) : m_lhs(lhs), m_rhs(rhs), m_functor(func) @@ -130,6 +132,10 @@ class CwiseBinaryOp : ei_no_assignment_operator, return m_functor.packetOp(m_lhs.template packet(index), m_rhs.template packet(index)); } + const _LhsNested& lhs() const { return m_lhs; } + const _RhsNested& rhs() const { return m_rhs; } + const BinaryOp& functor() const { return m_functor; } + protected: const LhsNested m_lhs; const RhsNested m_rhs; diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index a554de8e3..6ed96e286 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -457,6 +457,21 @@ template class MatrixBase void transposeInPlace(); const AdjointReturnType adjoint() const; void adjointInPlace(); + #ifndef EIGEN_NO_DEBUG + template + Derived& lazyAssign(const Transpose& other); + template + Derived& lazyAssign(const CwiseBinaryOp,Transpose,DerivedB>& other); + template + Derived& lazyAssign(const CwiseBinaryOp,DerivedA,Transpose >& other); + + template + Derived& lazyAssign(const CwiseUnaryOp, NestByValue > >& other); + template + Derived& lazyAssign(const CwiseBinaryOp,CwiseUnaryOp, NestByValue > >,DerivedB>& other); + template + Derived& lazyAssign(const CwiseBinaryOp,DerivedA,CwiseUnaryOp, NestByValue > > >& other); + #endif RowXpr row(int i); const RowXpr row(int i) const; diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 0787daa61..8527edc2a 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -267,4 +267,92 @@ inline void MatrixBase::adjointInPlace() derived() = adjoint().eval(); } +#ifndef EIGEN_NO_DEBUG + +// The following is to detect aliasing problems in the following common cases: +// a = a.transpose() +// a = a.transpose() + X +// a = X + a.transpose() +// a = a.adjoint() +// a = a.adjoint() + X +// a = X + a.adjoint() + +template::ActualAccess> +struct ei_extract_data_selector { + static typename T::Scalar* run(const T& m) + { + return &ei_blas_traits::extract(m).const_cast_derived().coeffRef(0,0); + } +}; + +template +struct ei_extract_data_selector { + static typename T::Scalar* run(const T&) { return 0; } +}; + +template typename T::Scalar* ei_extract_data(const T& m) +{ + return ei_extract_data_selector::run(m); +} + +template +template +Derived& MatrixBase::lazyAssign(const Transpose& other) +{ + ei_assert(ei_extract_data(other) != ei_extract_data(derived()) + && "aliasing detected during tranposition, please use transposeInPlace()"); + return lazyAssign(static_cast >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,Transpose,DerivedB>& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.lhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,Transpose,DerivedB> >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,DerivedA,Transpose >& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.rhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,DerivedA,Transpose > >& >(other)); +} + +template +template Derived& +MatrixBase:: +lazyAssign(const CwiseUnaryOp, NestByValue > >& other) +{ + ei_assert(ei_extract_data(other) != ei_extract_data(derived()) + && "aliasing detected during tranposition, please use adjointInPlace()"); + return lazyAssign(static_cast, NestByValue > > >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,CwiseUnaryOp, NestByValue > >,DerivedB>& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.lhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,CwiseUnaryOp, NestByValue > >,DerivedB> >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,DerivedA,CwiseUnaryOp, NestByValue > > >& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.rhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,DerivedA,CwiseUnaryOp, NestByValue > > > >& >(other)); +} +#endif + #endif // EIGEN_TRANSPOSE_H diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 300650338..964658c65 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -21,7 +21,7 @@ // You should have received a copy of the GNU Lesser General Public // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#define EIGEN_NO_ASSERTION_CHECKING + #include "main.h" template void adjoint(const MatrixType& m) @@ -110,7 +110,6 @@ template void adjoint(const MatrixType& m) m3.transposeInPlace(); VERIFY_IS_APPROX(m3,m1.conjugate()); - } void test_adjoint() @@ -125,5 +124,16 @@ void test_adjoint() } // test a large matrix only once CALL_SUBTEST( adjoint(Matrix()) ); + + { + MatrixXcf a(10,10), b(10,10); + VERIFY_RAISES_ASSERT(a = a.transpose()); + VERIFY_RAISES_ASSERT(a = a.transpose() + b); + VERIFY_RAISES_ASSERT(a = b + a.transpose()); + VERIFY_RAISES_ASSERT(a = a.conjugate().transpose()); + VERIFY_RAISES_ASSERT(a = a.adjoint()); + VERIFY_RAISES_ASSERT(a = a.adjoint() + b); + VERIFY_RAISES_ASSERT(a = b + a.adjoint()); + } } From 03c1e79f352ee3a815c2b9f0e07842cd2500d017 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sat, 15 Aug 2009 19:20:48 -0400 Subject: [PATCH 45/58] svd: sort in decreasing order, remove unused code --- Eigen/src/SVD/SVD.h | 42 +----------------------------------------- test/svd.cpp | 1 - 2 files changed, 1 insertion(+), 42 deletions(-) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index acc567f94..ddc0de4ce 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -98,7 +98,6 @@ template class SVD } void compute(const MatrixType& matrix); - SVD& sort(); template void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const; @@ -392,7 +391,7 @@ void SVD::compute(const MatrixType& matrix) for (int i=0; i::compute(const MatrixType& matrix) m_isInitialized = true; } -template -SVD& SVD::sort() -{ - ei_assert(m_isInitialized && "SVD is not initialized."); - - int mu = m_matU.rows(); - int mv = m_matV.rows(); - int n = m_matU.cols(); - - for (int i=0; i p) - { - k = j; - p = m_sigma.coeff(j); - } - } - if (k != i) - { - m_sigma.coeffRef(k) = m_sigma.coeff(i); // i.e. - m_sigma.coeffRef(i) = p; // swaps the i-th and the k-th elements - - int j = mu; - for(int s=0; j!=0; ++s, --j) - std::swap(m_matU.coeffRef(s,i), m_matU.coeffRef(s,k)); - - j = mv; - for (int s=0; j!=0; ++s, --j) - std::swap(m_matV.coeffRef(s,i), m_matV.coeffRef(s,k)); - } - } - return *this; -} - /** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A. * The parts of the solution corresponding to zero singular values are ignored. * diff --git a/test/svd.cpp b/test/svd.cpp index 93e5ea017..2ccd94764 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -95,7 +95,6 @@ template void svd_verify_assert() VERIFY_RAISES_ASSERT(svd.matrixU()) VERIFY_RAISES_ASSERT(svd.singularValues()) VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.sort()) VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp)) VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp)) VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp)) From 044dd0c1ddaee5d432c81b8d5d8c3674b7a003f2 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 23:37:20 +0200 Subject: [PATCH 46/58] revert previous change in Quaternion::setFromTwoVectors --- Eigen/Geometry | 1 + Eigen/src/Geometry/Quaternion.h | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Eigen/Geometry b/Eigen/Geometry index 9cae3459c..c931f28fe 100644 --- a/Eigen/Geometry +++ b/Eigen/Geometry @@ -7,6 +7,7 @@ #include "Array" #include "SVD" +#include "LU" #include #ifndef M_PI diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 3d59533f4..2f9f97807 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -373,7 +373,7 @@ inline Quaternion& Quaternion::setFromTwoVectors(const MatrixBas c = std::max(c,-1); Matrix m; m << v0.transpose(), v1.transpose(); SVD > svd(m); - Vector3 axis = svd.matrixV().col(0); + Vector3 axis = svd.matrixV().col(2); Scalar w2 = (Scalar(1)+c)*Scalar(0.5); this->w() = ei_sqrt(w2); From f5f2b222a3862d5ccc19e6e44e4f37de1a56dfb5 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 16 Aug 2009 00:02:36 +0200 Subject: [PATCH 47/58] make SVD reuses applyJacobi --- Eigen/src/SVD/SVD.h | 26 +++++--------------------- 1 file changed, 5 insertions(+), 21 deletions(-) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index ddc0de4ce..2948087a2 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -307,13 +307,7 @@ void SVD::compute(const MatrixType& matrix) h = Scalar(1.0)/h; c = g*h; s = -f*h; - for (j=0; j::compute(const MatrixType& matrix) y = W[i]; h = s*g; g = c*g; + z = pythag(f,h); rv1[j] = z; c = f/z; @@ -354,13 +349,8 @@ void SVD::compute(const MatrixType& matrix) g = g*c - x*s; h = y*s; y *= c; - for (jj=0; jj::compute(const MatrixType& matrix) } f = c*g + s*y; x = c*y - s*g; - for (jj=0; jj Date: Sun, 16 Aug 2009 00:14:05 +0200 Subject: [PATCH 48/58] rename back MayAliasBit to EvalBeforeAssigningBit --- Eigen/src/Core/Assign.h | 2 +- Eigen/src/Core/Flagged.h | 4 ++-- Eigen/src/Core/MapBase.h | 4 ++-- Eigen/src/Core/Matrix.h | 2 +- Eigen/src/Core/MatrixBase.h | 8 ++++---- Eigen/src/Core/ProductBase.h | 8 ++++---- Eigen/src/Core/TriangularMatrix.h | 6 +++--- Eigen/src/Core/products/GeneralUnrolled.h | 2 +- Eigen/src/Core/util/Constants.h | 8 +++----- 9 files changed, 21 insertions(+), 23 deletions(-) diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index b0e69224e..4bd1046a7 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -434,7 +434,7 @@ EIGEN_STRONG_INLINE Derived& MatrixBase } template::marked() const /** \deprecated use MatrixBase::noalias() * - * \returns an expression of *this with the MayAliasBit flag removed. + * \returns an expression of *this with the EvalBeforeAssigningBit flag removed. * * Example: \include MatrixBase_lazy.cpp * Output: \verbinclude MatrixBase_lazy.out @@ -140,7 +140,7 @@ MatrixBase::marked() const * \sa class Flagged, marked() */ template -inline const Flagged +inline const Flagged MatrixBase::lazy() const { return derived(); diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index b9e9374be..9cb085b61 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -179,11 +179,11 @@ template class MapBase // explicitly add these two overloads. // Maybe there exists a better solution though. template - Derived& operator+=(const Flagged, 0, MayAliasBit>& other) + Derived& operator+=(const Flagged, 0, EvalBeforeAssigningBit>& other) { return Base::operator+=(other); } template - Derived& operator-=(const Flagged, 0, MayAliasBit>& other) + Derived& operator-=(const Flagged, 0, EvalBeforeAssigningBit>& other) { return Base::operator-=(other); } template diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index d129535b6..f58424ba2 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -577,7 +577,7 @@ class Matrix template EIGEN_STRONG_INLINE Matrix& _set(const MatrixBase& other) { - _set_selector(other.derived(), typename ei_meta_if(int(OtherDerived::Flags) & MayAliasBit), ei_meta_true, ei_meta_false>::ret()); + _set_selector(other.derived(), typename ei_meta_if(int(OtherDerived::Flags) & EvalBeforeAssigningBit), ei_meta_true, ei_meta_false>::ret()); return *this; } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 6ed96e286..380ac89c0 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -329,7 +329,7 @@ template class MatrixBase /** \deprecated because .lazy() is deprecated * Overloaded for cache friendly product evaluation */ template - Derived& lazyAssign(const Flagged& other) + Derived& lazyAssign(const Flagged& other) { return lazyAssign(other._expression()); } template @@ -337,11 +337,11 @@ template class MatrixBase template Derived& operator+=(const Flagged, 0, - MayAliasBit>& other); + EvalBeforeAssigningBit>& other); template Derived& operator-=(const Flagged, 0, - MayAliasBit>& other); + EvalBeforeAssigningBit>& other); #endif // not EIGEN_PARSED_BY_DOXYGEN CommaInitializer operator<< (const Scalar& s); @@ -630,7 +630,7 @@ template class MatrixBase template const Flagged marked() const; - const Flagged lazy() const; + const Flagged lazy() const; NoAlias noalias(); diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index 867aa7a15..5e259b7f7 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -39,7 +39,7 @@ struct ei_traits > ColsAtCompileTime = ei_traits::ColsAtCompileTime, MaxRowsAtCompileTime = ei_traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = ei_traits::MaxColsAtCompileTime, - Flags = EvalBeforeNestingBit | MayAliasBit, + Flags = EvalBeforeNestingBit | EvalBeforeAssigningBit, CoeffReadCost = 0 // FIXME why is it needed ? }; }; @@ -119,7 +119,7 @@ class ProductBase : public MatrixBase return res; } - const Flagged lazy() const + const Flagged lazy() const { return *this; } @@ -228,7 +228,7 @@ Derived& MatrixBase::lazyAssign(const ProductBase template Derived& MatrixBase::operator+=(const Flagged, 0, - MayAliasBit>& other) + EvalBeforeAssigningBit>& other) { other._expression().derived().addTo(derived()); return derived(); } @@ -238,7 +238,7 @@ Derived& MatrixBase::operator+=(const Flagged template Derived& MatrixBase::operator-=(const Flagged, 0, - MayAliasBit>& other) + EvalBeforeAssigningBit>& other) { other._expression().derived().subTo(derived()); return derived(); } diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index cfaf89733..b0362f20c 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -477,7 +477,7 @@ template inline TriangularView& TriangularView::operator=(const MatrixBase& other) { - if(OtherDerived::Flags & MayAliasBit) + if(OtherDerived::Flags & EvalBeforeAssigningBit) { typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); other_evaluated.template triangularView().lazyAssign(other.derived()); @@ -512,7 +512,7 @@ inline TriangularView& TriangularView::operator=(const TriangularBase& other) { ei_assert(Mode == OtherDerived::Mode); - if(ei_traits::Flags & MayAliasBit) + if(ei_traits::Flags & EvalBeforeAssigningBit) { typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); other_evaluated.template triangularView().lazyAssign(other.derived()); @@ -548,7 +548,7 @@ template template void TriangularBase::evalToDense(MatrixBase &other) const { - if(ei_traits::Flags & MayAliasBit) + if(ei_traits::Flags & EvalBeforeAssigningBit) { typename Derived::PlainMatrixType other_evaluated(rows(), cols()); evalToDenseLazy(other_evaluated); diff --git a/Eigen/src/Core/products/GeneralUnrolled.h b/Eigen/src/Core/products/GeneralUnrolled.h index 2961d9452..7241976a8 100644 --- a/Eigen/src/Core/products/GeneralUnrolled.h +++ b/Eigen/src/Core/products/GeneralUnrolled.h @@ -76,7 +76,7 @@ struct ei_traits > RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) - | MayAliasBit + | EvalBeforeAssigningBit | EvalBeforeNestingBit | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0) | (LhsFlags & RhsFlags & AlignedBit), diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 1d55217f5..21e7f62c3 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -76,10 +76,8 @@ const unsigned int EvalBeforeNestingBit = 0x2; /** \ingroup flags * - * Means the expression cannot be evaluated safely if the result alias one - * of the operands of the expression. Therefore it should be evaluated - * before any assignement. */ -const unsigned int MayAliasBit = 0x4; + * means the expression should be evaluated before any assignement */ +const unsigned int EvalBeforeAssigningBit = 0x4; /** \ingroup flags * @@ -184,7 +182,7 @@ const unsigned int SparseBit = 0x1000; // list of flags that are inherited by default const unsigned int HereditaryBits = RowMajorBit | EvalBeforeNestingBit - | MayAliasBit + | EvalBeforeAssigningBit | SparseBit; // Possible values for the Mode parameter of part() From ee982709d37e219e9d67beaf777a5bf2131e835e Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sat, 15 Aug 2009 23:12:39 -0400 Subject: [PATCH 49/58] in all decs, make the compute() methods return *this (implements feature request #18) --- Eigen/src/Cholesky/LDLT.h | 7 ++++--- Eigen/src/Cholesky/LLT.h | 8 ++++++-- Eigen/src/LU/LU.h | 7 +++++-- Eigen/src/LU/PartialLU.h | 5 +++-- Eigen/src/QR/EigenSolver.h | 5 +++-- Eigen/src/QR/QR.h | 5 +++-- Eigen/src/QR/SelfAdjointEigenSolver.h | 12 +++++++----- Eigen/src/SVD/JacobiSquareSVD.h | 5 +++-- Eigen/src/SVD/SVD.h | 9 ++++++--- 9 files changed, 40 insertions(+), 23 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 1cb1294ac..adca9fe2e 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -123,7 +123,7 @@ template class LDLT template bool solveInPlace(MatrixBase &bAndX) const; - void compute(const MatrixType& matrix); + LDLT& compute(const MatrixType& matrix); protected: /** \internal @@ -142,7 +142,7 @@ template class LDLT /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix */ template -void LDLT::compute(const MatrixType& a) +LDLT& LDLT::compute(const MatrixType& a) { ei_assert(a.rows()==a.cols()); const int size = a.rows(); @@ -154,7 +154,7 @@ void LDLT::compute(const MatrixType& a) m_transpositions.setZero(); m_sign = ei_real(a.coeff(0,0))>0 ? 1:-1; m_isInitialized = true; - return; + return *this; } RealScalar cutoff = 0, biggest_in_corner; @@ -235,6 +235,7 @@ void LDLT::compute(const MatrixType& a) } m_isInitialized = true; + return *this; } /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index ef04b8fe4..3ce85b2d9 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -105,7 +105,7 @@ template class LLT template bool solveInPlace(MatrixBase &bAndX) const; - void compute(const MatrixType& matrix); + LLT& compute(const MatrixType& matrix); protected: /** \internal @@ -213,9 +213,12 @@ template struct LLT_Traits }; /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix + * + * + * \returns a reference to *this */ template -void LLT::compute(const MatrixType& a) +LLT& LLT::compute(const MatrixType& a) { assert(a.rows()==a.cols()); const int size = a.rows(); @@ -223,6 +226,7 @@ void LLT::compute(const MatrixType& a) m_matrix = a; m_isInitialized = Traits::inplace_decomposition(m_matrix); + return *this; } /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index 733fa0cbc..f072ff906 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -111,8 +111,10 @@ template class LU * * \param matrix the matrix of which to compute the LU decomposition. * It is required to be nonzero. + * + * \returns a reference to *this */ - void compute(const MatrixType& matrix); + LU& compute(const MatrixType& matrix); /** \returns the LU decomposition matrix: the upper-triangular part is U, the * unit-lower-triangular part is L (at least for square matrices; in the non-square @@ -377,7 +379,7 @@ LU::LU(const MatrixType& matrix) } template -void LU::compute(const MatrixType& matrix) +LU& LU::compute(const MatrixType& matrix) { m_originalMatrix = &matrix; m_lu = matrix; @@ -447,6 +449,7 @@ void LU::compute(const MatrixType& matrix) std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k))); m_det_pq = (number_of_transpositions%2) ? -1 : 1; + return *this; } template diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h index f5c9b8ae9..eb5c3d34b 100644 --- a/Eigen/src/LU/PartialLU.h +++ b/Eigen/src/LU/PartialLU.h @@ -87,7 +87,7 @@ template class PartialLU */ PartialLU(const MatrixType& matrix); - void compute(const MatrixType& matrix); + PartialLU& compute(const MatrixType& matrix); /** \returns the LU decomposition matrix: the upper-triangular part is U, the * unit-lower-triangular part is L (at least for square matrices; in the non-square @@ -350,7 +350,7 @@ void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, int& n } template -void PartialLU::compute(const MatrixType& matrix) +PartialLU& PartialLU::compute(const MatrixType& matrix) { m_lu = matrix; m_p.resize(matrix.rows()); @@ -369,6 +369,7 @@ void PartialLU::compute(const MatrixType& matrix) std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); m_isInitialized = true; + return *this; } template diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h index dbe1aca89..bd7a76c49 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/Eigen/src/QR/EigenSolver.h @@ -118,7 +118,7 @@ template class EigenSolver return m_eivalues; } - void compute(const MatrixType& matrix); + EigenSolver& compute(const MatrixType& matrix); private: @@ -189,7 +189,7 @@ typename EigenSolver::EigenvectorType EigenSolver::eigen } template -void EigenSolver::compute(const MatrixType& matrix) +EigenSolver& EigenSolver::compute(const MatrixType& matrix) { assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); @@ -205,6 +205,7 @@ void EigenSolver::compute(const MatrixType& matrix) hqr2(matH); m_isInitialized = true; + return *this; } // Nonsymmetric reduction to Hessenberg form. diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 98ef4da25..d37a019ff 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -101,7 +101,7 @@ template class HouseholderQR */ const MatrixType& matrixQR() const { return m_qr; } - void compute(const MatrixType& matrix); + HouseholderQR& compute(const MatrixType& matrix); protected: MatrixType m_qr; @@ -112,7 +112,7 @@ template class HouseholderQR #ifndef EIGEN_HIDE_HEAVY_CODE template -void HouseholderQR::compute(const MatrixType& matrix) +HouseholderQR& HouseholderQR::compute(const MatrixType& matrix) { m_qr = matrix; m_hCoeffs.resize(matrix.cols()); @@ -175,6 +175,7 @@ void HouseholderQR::compute(const MatrixType& matrix) } } m_isInitialized = true; + return *this; } template diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index d07d28424..809a717b9 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -90,9 +90,9 @@ template class SelfAdjointEigenSolver compute(matA, matB, computeEigenvectors); } - void compute(const MatrixType& matrix, bool computeEigenvectors = true); + SelfAdjointEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); - void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true); + SelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true); /** \returns the computed eigen vectors as a matrix of column vectors */ MatrixType eigenvectors(void) const @@ -182,7 +182,7 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool) */ template -void SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool computeEigenvectors) +SelfAdjointEigenSolver& SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool computeEigenvectors) { #ifndef NDEBUG m_eigenvectorsOk = computeEigenvectors; @@ -195,7 +195,7 @@ void SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool { m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0)); m_eivec.setOnes(); - return; + return *this; } m_eivec = matrix; @@ -240,6 +240,7 @@ void SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool m_eivec.col(i).swap(m_eivec.col(k+i)); } } + return *this; } /** Computes the eigenvalues of the generalized eigen problem @@ -250,7 +251,7 @@ void SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool) */ template -void SelfAdjointEigenSolver:: +SelfAdjointEigenSolver& SelfAdjointEigenSolver:: compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors) { ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); @@ -282,6 +283,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors for (int i=0; i class JacobiSquareSV compute(matrix); } - void compute(const MatrixType& matrix); + JacobiSquareSVD& compute(const MatrixType& matrix); const MatrixUType& matrixU() const { @@ -95,7 +95,7 @@ template class JacobiSquareSV }; template -void JacobiSquareSVD::compute(const MatrixType& matrix) +JacobiSquareSVD& JacobiSquareSVD::compute(const MatrixType& matrix) { MatrixType work_matrix(matrix); int size = matrix.rows(); @@ -164,5 +164,6 @@ sweep_again: } m_isInitialized = true; + return *this; } #endif // EIGEN_JACOBISQUARESVD_H diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 2948087a2..1a7e6c49a 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -97,7 +97,7 @@ template class SVD return m_matV; } - void compute(const MatrixType& matrix); + SVD& compute(const MatrixType& matrix); template void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const; @@ -138,9 +138,11 @@ template class SVD /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix * * \note this code has been adapted from Numerical Recipes, third edition. + * + * \returns a reference to *this */ template -void SVD::compute(const MatrixType& matrix) +SVD& SVD::compute(const MatrixType& matrix) { const int m = matrix.rows(); const int n = matrix.cols(); @@ -157,7 +159,7 @@ void SVD::compute(const MatrixType& matrix) SingularValuesType& W = m_sigma; bool flag; - int i,its,j,jj,k,l,nm; + int i,its,j,k,l,nm; Scalar anorm, c, f, g, h, s, scale, x, y, z; bool convergence = true; Scalar eps = precision(); @@ -392,6 +394,7 @@ void SVD::compute(const MatrixType& matrix) m_matU = A.block(0,0,m,m); m_isInitialized = true; + return *this; } /** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A. From fc9480cbb38db2c7585e3ac5d14719237a85d2d7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 16 Aug 2009 10:55:10 +0200 Subject: [PATCH 50/58] bugfix in compute_matrix_flags, optimization in LU, improve doc, and workaround aliasing detection in MatrixBase_eval snippet (not very nice but I don't know how to do it in a better way) --- Eigen/src/Core/NoAlias.h | 52 +++++++---- Eigen/src/Core/ProductBase.h | 6 +- Eigen/src/Core/util/XprHelper.h | 2 +- Eigen/src/LU/LU.h | 2 +- Eigen/src/LU/PartialLU.h | 2 +- doc/I01_TopicLazyEvaluation.dox | 2 +- doc/I02_HiPerformance.dox | 154 +++++++------------------------ doc/snippets/MatrixBase_eval.cpp | 2 +- 8 files changed, 77 insertions(+), 145 deletions(-) diff --git a/Eigen/src/Core/NoAlias.h b/Eigen/src/Core/NoAlias.h index 5c2b6e3ca..a45493ddc 100644 --- a/Eigen/src/Core/NoAlias.h +++ b/Eigen/src/Core/NoAlias.h @@ -31,8 +31,9 @@ * * \param ExpressionType the type of the object on which to do the lazy assignment * - * This class represents an expression with a special assignment operator (operator=) + * This class represents an expression with special assignment operators * assuming no aliasing between the target expression and the source expression. + * More precisely it alloas to bypass the EvalBeforeAssignBit flag of the source expression. * It is the return type of MatrixBase::noalias() * and most of the time this is the only way it is used. * @@ -44,44 +45,61 @@ class NoAlias public: NoAlias(ExpressionType& expression) : m_expression(expression) {} - /** Behaves like MatrixBase::lazyAssign() */ + /** Behaves like MatrixBase::lazyAssign(other) + * \sa MatrixBase::lazyAssign() */ template ExpressionType& operator=(const MatrixBase& other) - { - return m_expression.lazyAssign(other.derived()); - } + { return m_expression.lazyAssign(other.derived()); } - // TODO could be removed if we decide that += is noalias by default + /** \sa MatrixBase::operator+= */ template ExpressionType& operator+=(const MatrixBase& other) - { - return m_expression.lazyAssign(m_expression + other.derived()); - } + { return m_expression.lazyAssign(m_expression + other.derived()); } - // TODO could be removed if we decide that += is noalias by default + /** \sa MatrixBase::operator-= */ template ExpressionType& operator-=(const MatrixBase& other) - { - return m_expression.lazyAssign(m_expression - other.derived()); - } + { return m_expression.lazyAssign(m_expression - other.derived()); } - // TODO could be removed if we decide that += is noalias by default +#ifndef EIGEN_PARSED_BY_DOXYGEN template ExpressionType& operator+=(const ProductBase& other) { other.derived().addTo(m_expression); return m_expression; } - // TODO could be removed if we decide that += is noalias by default template ExpressionType& operator-=(const ProductBase& other) { other.derived().subTo(m_expression); return m_expression; } +#endif protected: ExpressionType& m_expression; }; - /** \returns a pseudo expression of \c *this with an operator= assuming - * no aliasing between \c *this and the source expression + * no aliasing between \c *this and the source expression. + * + * More precisely, noalias() allows to bypass the EvalBeforeAssignBit flag. + * Currently, even though several expressions may alias, only product + * expressions have this flag. Therefore, noalias() is only usefull when + * the source expression contains a matrix product. + * + * Here are some examples where noalias is usefull: + * \code + * D.noalias() = A * B; + * D.noalias() += A.transpose() * B; + * D.noalias() -= 2 * A * B.adjoint(); + * \endcode + * + * On the other hand the following example will lead to a \b wrong result: + * \code + * A.noalias() = A * B; + * \endcode + * because the result matrix A is also an operand of the matrix product. Therefore, + * there is no alternative than evaluating A * B in a temporary, that is the default + * behavior when you write: + * \code + * A = A * B; + * \endcode * * \sa class NoAlias */ diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index 5e259b7f7..b2c4cd989 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -119,10 +119,8 @@ class ProductBase : public MatrixBase return res; } - const Flagged lazy() const - { - return *this; - } + EIGEN_DEPRECATED const Flagged lazy() const + { return *this; } const _LhsNested& lhs() const { return m_lhs; } const _RhsNested& rhs() const { return m_rhs; } diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index d392e27ba..871259b08 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -90,7 +90,7 @@ class ei_compute_matrix_flags : MaxRows==1 ? MaxCols : row_major_bit ? MaxCols : MaxRows, is_big = inner_max_size == Dynamic, - is_packet_size_multiple = Rows==Dynamic || Cols==Dynamic || ((Cols*Rows) % ei_packet_traits::size) == 0, + is_packet_size_multiple = MaxRows==Dynamic || MaxCols==Dynamic || ((MaxCols*MaxRows) % ei_packet_traits::size) == 0, aligned_bit = (((Options&DontAlign)==0) && (is_big || is_packet_size_multiple)) ? AlignedBit : 0, packet_access_bit = ei_packet_traits::size > 1 && aligned_bit ? PacketAccessBit : 0 }; diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index f072ff906..d32fe08a1 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -437,7 +437,7 @@ LU& LU::compute(const MatrixType& matrix) if(k \code m1 += m2 * m3; \endcode \code temp = m2 * m3; m1 += temp; \endcode -\code m1 += (m2 * m3).lazy(); \endcode -Use .lazy() to tell Eigen the result and right-hand-sides do not alias. +\code m1.noalias() += m2 * m3; \endcode +Use .noalias() to tell Eigen the result and right-hand-sides do not alias. + Otherwise the product m2 * m3 is evaluated into a temporary. -\code m1 += (s1 * (m2 * m3)).lazy(); \endcode -\code temp = (m2 * m3).lazy(); m1 += s1 * temp; \endcode -\code m1 += (s1 * m2 * m3).lazy(); \endcode -This is because m2 * m3 is immediately evaluated by the scalar product.
- Make sure the matrix product is the top most expression. + + +\code m1.noalias() += s1 * (m2 * m3); \endcode +This is a special feature of Eigen. Here the product between a scalar + and a matrix product does not evaluate the matrix product but instead it + returns a matrix product expression tracking the scalar scaling factor.
+ Without this optimization, the matrix product would be evaluated into a + temporary as in the next example. -\code m1 += s1 * (m2 * m3).lazy(); \endcode -\code m1 += s1 * m2 * m3; // using a naive product \endcode -\code m1 += (s1 * m2 * m3).lazy(); \endcode -Even though this expression is evaluated without temporary, it is actually even - worse than the previous case because here the .lazy() enforces Eigen to use a - naive (and slow) evaluation of the product. +\code m1.noalias() += (m2 * m3).transpose(); \endcode +\code temp = m2 * m3; m1 += temp.transpose(); \endcode +\code m1.noalias() += m3.adjoint() * m3.adjoint(); \endcode +This is because the product expression has the EvalBeforeNesting bit which + enforces the evaluation of the product by the Tranpose expression. \code m1 = m1 + m2 * m3; \endcode @@ -78,112 +81,25 @@ handled by a single GEMM-like call are correctly detected. and so the matrix product will be immediately evaluated. -\code m1 += ((s1*m2).block(....) * m3).lazy(); \endcode -\code temp = (s1*m2).block(....); m1 += (temp * m3).lazy(); \endcode -\code m1 += (s1 * m2.block(....) * m3).lazy(); \endcode +\code m1.noalias() = m4 + m2 * m3; \endcode +\code temp = m2 * m3; m1 = m4 + temp; \endcode +\code m1 = m4; m1.noalias() += m2 * m3; \endcode +First of all, here the .noalias() in the first expression is useless because + m2*m3 will be evaluated anyway. However, note how this expression can be rewritten + so that no temporary is evaluated. (tip: for very small fixed size matrix + it is slighlty better to rewrite it like this: m1.noalias() = m2 * m3; m1 += m4; + + +\code m1.noalias() += ((s1*m2).block(....) * m3); \endcode +\code temp = (s1*m2).block(....); m1 += temp * m3; \endcode +\code m1.noalias() += s1 * m2.block(....) * m3; \endcode This is because our expression analyzer is currently not able to extract trivial expressions nested in a Block expression. Therefore the nested scalar multiple cannot be properly extracted. -Of course all these remarks hold for all other kind of products that we will describe in the following paragraphs. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
BLAS equivalent routineEfficient version
(compile to a single optimized evaluation)		Less efficient equivalent version
(requires multiple evaluations)		comments
GEMMm1 = s1 * m2 * m3m1 = s1 * (m2 * m3)This is because m2 * m3 is evaluated by the scalar product.
GEMMm1 += s1 * m2.adjoint() * m3m1 += (s1 * m2).adjoint() * m3This is because our expression analyzer stops at the first transpose expression and cannot extract the nested scalar multiple.
GEMMm1 += m2.adjoint() * m3m1 += m2.conjugate().transpose() * m3For the same reason. Use .adjoint() or .transpose().conjugate()
GEMMm1 -= (-(s0*m2).conjugate()*s1) * (s2 * m3.adjoint() * s3)Note that s0 is automatically conjugated during the simplification of the expression.
SYRm.sefadjointView().rankUpdate(v,s)Computes m += s * v * v.adjoint()
SYR2m.sefadjointView().rankUpdate(u,v,s)Computes m += s * u * v.adjoint() + s * v * u.adjoint()
SYRKm1.sefadjointView().rankUpdate(m2.adjoint(),s)Computes m1 += s * m2.adjoint() * m2
SYMM/HEMMm3 -= s1 * m1.sefadjointView() * m2.adjoint()
SYMM/HEMMm3 += s1 * m2.transpose() * m1.conjugate().sefadjointView()
TRMMm3 -= s1 * m1.triangularView() * m2.adjoint()
TRSV / TRSMm1.adjoint().triangularView().solveInPlace(m2)
+Of course all these remarks hold for all other kind of products involving triangular or selfadjoint matrices. */ diff --git a/doc/snippets/MatrixBase_eval.cpp b/doc/snippets/MatrixBase_eval.cpp index 732293043..d70424562 100644 --- a/doc/snippets/MatrixBase_eval.cpp +++ b/doc/snippets/MatrixBase_eval.cpp @@ -4,7 +4,7 @@ m = M; cout << "Here is the matrix m:" << endl << m << endl; cout << "Now we want to replace m by its own transpose." << endl; cout << "If we do m = m.transpose(), then m becomes:" << endl; -m = m.transpose(); +m = m.transpose() * 1; cout << m << endl << "which is wrong!" << endl; cout << "Now let us instead do m = m.transpose().eval(). Then m becomes" << endl; m = M; From 5274c5c3262cea98d5b78962af5680e4ae6050c0 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 16 Aug 2009 11:01:32 +0200 Subject: [PATCH 51/58] quick update of TopicLazyEvaluation --- doc/I01_TopicLazyEvaluation.dox | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/doc/I01_TopicLazyEvaluation.dox b/doc/I01_TopicLazyEvaluation.dox index 3e2b9db5f..f40afa06d 100644 --- a/doc/I01_TopicLazyEvaluation.dox +++ b/doc/I01_TopicLazyEvaluation.dox @@ -42,11 +42,11 @@ Eigen chooses lazy evaluation at every stage in that example, which is clearly t Eigen first evaluates matrix * matrix into a temporary matrix, and then copies it into the original \c matrix. This guarantees a correct result as we saw above that lazy evaluation gives wrong results with matrix products. It also doesn't cost much, as the cost of the matrix product itself is much higher. -What if you know what you are doing and want to force lazy evaluation? Then use \link MatrixBase::lazy() .lazy()\endlink instead. Here is an example: +What if you know that the result does no alias the operand of the product and want to force lazy evaluation? Then use \link MatrixBase::noalias() .noalias()\endlink instead. Here is an example: -\code matrix1 = (matrix2 * matrix2).lazy(); \endcode +\code matrix1.noalias() = matrix2 * matrix2; \endcode -Here, since we know that matrix2 is not the same matrix as matrix1, we know that lazy evaluation is not dangerous, so we may force lazy evaluation. Concretely, the effect of lazy() here is to remove the evaluate-before-assigning \link flags flag\endlink and also the evaluate-before-nesting \link flags flag\endlink which we now discuss. +Here, since we know that matrix2 is not the same matrix as matrix1, we know that lazy evaluation is not dangerous, so we may force lazy evaluation. Concretely, the effect of noalias() here is to bypass the evaluate-before-assigning \link flags flag\endlink. The second circumstance in which Eigen chooses immediate evaluation, is when it sees a nested expression such as a + b where \c b is already an expression having the evaluate-before-nesting \link flags flag\endlink. Again, the most important example of such an expression is the \link Product matrix product expression\endlink. For example, when you do @@ -54,8 +54,6 @@ Here, since we know that matrix2 is not the same matrix as matrix1, we know that the product matrix3 * matrix4 gets evaluated immediately into a temporary matrix. Indeed, experiments showed that it is often beneficial for performance to evaluate immediately matrix products when they are nested into bigger expressions. -Again, \link MatrixBase::lazy() .lazy()\endlink can be used to force lazy evaluation here. - The third circumstance in which Eigen chooses immediate evaluation, is when its cost model shows that the total cost of an operation is reduced if a sub-expression gets evaluated into a temporary. Indeed, in certain cases, an intermediate result is sufficiently costly to compute and is reused sufficiently many times, that is worth "caching". Here is an example: \code matrix1 = matrix2 * (matrix3 + matrix4); \endcode From ef13c3e754008e507eb67a9c5eab8dcb812777a1 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 16 Aug 2009 11:51:46 +0200 Subject: [PATCH 52/58] add normalize and normalized overloads in AlignedVector3 --- unsupported/Eigen/AlignedVector3 | 11 +++++++++++ unsupported/test/alignedvector3.cpp | 6 +++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/unsupported/Eigen/AlignedVector3 b/unsupported/Eigen/AlignedVector3 index f9f364eba..854f0601f 100644 --- a/unsupported/Eigen/AlignedVector3 +++ b/unsupported/Eigen/AlignedVector3 @@ -152,8 +152,19 @@ template class AlignedVector3 { ei_assert(m_coeffs.w()==Scalar(0)); ei_assert(other.m_coeffs.w()==Scalar(0)); + Scalar r = m_coeffs.dot(other.m_coeffs); return m_coeffs.dot(other.m_coeffs); } + + inline void normalize() + { + m_coeffs /= norm(); + } + + inline AlignedVector3 normalized() + { + return AlignedVector3(m_coeffs / norm()); + } inline Scalar sum() const { diff --git a/unsupported/test/alignedvector3.cpp b/unsupported/test/alignedvector3.cpp index cf33dcceb..52d39ce50 100644 --- a/unsupported/test/alignedvector3.cpp +++ b/unsupported/test/alignedvector3.cpp @@ -40,7 +40,7 @@ void alignedvector3() VERIFY_IS_APPROX(f1,r1); VERIFY_IS_APPROX(f4,r4); - + VERIFY_IS_APPROX(f4+f1,r4+r1); VERIFY_IS_APPROX(f4-f1,r4-r1); VERIFY_IS_APPROX(f4+f1-f2,r4+r1-r2); @@ -56,6 +56,10 @@ void alignedvector3() VERIFY_IS_APPROX(f2.dot(f3),r2.dot(r3)); VERIFY_IS_APPROX(f2.cross(f3),r2.cross(r3)); VERIFY_IS_APPROX(f2.norm(),r2.norm()); + + VERIFY_IS_APPROX(f2.normalized(),r2.normalized()); + + VERIFY_IS_APPROX((f2+f1).normalized(),(r2+r1).normalized()); f2.normalize(); r2.normalize(); From 737bed19c1fdb01568706bca19666531dda681a7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 16 Aug 2009 19:22:15 +0200 Subject: [PATCH 53/58] make HouseholderQR uses the Householder module --- Eigen/QR | 1 + Eigen/src/Core/MatrixBase.h | 9 ++- Eigen/src/Core/NoAlias.h | 10 +-- Eigen/src/Core/Product.h | 4 +- Eigen/src/Householder/Householder.h | 50 ++++++++++---- Eigen/src/QR/QR.h | 103 ++++++++++------------------ test/qr.cpp | 21 +++--- 7 files changed, 96 insertions(+), 102 deletions(-) diff --git a/Eigen/QR b/Eigen/QR index ecebd32e4..151c0b31b 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -7,6 +7,7 @@ #include "Cholesky" #include "Jacobi" +#include "Householder" // Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module #if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 380ac89c0..9e92c043f 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -786,15 +786,18 @@ template class MatrixBase ////////// Householder module /////////// + void makeHouseholderInPlace(RealScalar *tau, Scalar *beta); template void makeHouseholder(EssentialPart *essential, - RealScalar *beta) const; + RealScalar *tau, Scalar *beta) const; template void applyHouseholderOnTheLeft(const EssentialPart& essential, - const RealScalar& beta); + const RealScalar& tau, + Scalar* workspace); template void applyHouseholderOnTheRight(const EssentialPart& essential, - const RealScalar& beta); + const RealScalar& tau, + Scalar* workspace); ///////// Jacobi module ///////// diff --git a/Eigen/src/Core/NoAlias.h b/Eigen/src/Core/NoAlias.h index a45493ddc..66d8d834d 100644 --- a/Eigen/src/Core/NoAlias.h +++ b/Eigen/src/Core/NoAlias.h @@ -48,26 +48,26 @@ class NoAlias /** Behaves like MatrixBase::lazyAssign(other) * \sa MatrixBase::lazyAssign() */ template - ExpressionType& operator=(const MatrixBase& other) + EIGEN_STRONG_INLINE ExpressionType& operator=(const MatrixBase& other) { return m_expression.lazyAssign(other.derived()); } /** \sa MatrixBase::operator+= */ template - ExpressionType& operator+=(const MatrixBase& other) + EIGEN_STRONG_INLINE ExpressionType& operator+=(const MatrixBase& other) { return m_expression.lazyAssign(m_expression + other.derived()); } /** \sa MatrixBase::operator-= */ template - ExpressionType& operator-=(const MatrixBase& other) + EIGEN_STRONG_INLINE ExpressionType& operator-=(const MatrixBase& other) { return m_expression.lazyAssign(m_expression - other.derived()); } #ifndef EIGEN_PARSED_BY_DOXYGEN template - ExpressionType& operator+=(const ProductBase& other) + EIGEN_STRONG_INLINE ExpressionType& operator+=(const ProductBase& other) { other.derived().addTo(m_expression); return m_expression; } template - ExpressionType& operator-=(const ProductBase& other) + EIGEN_STRONG_INLINE ExpressionType& operator-=(const ProductBase& other) { other.derived().subTo(m_expression); return m_expression; } #endif diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index fba4241b1..e1e106b80 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -187,7 +187,7 @@ class GeneralProduct template<> struct ei_outer_product_selector { template - static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) { + EIGEN_DONT_INLINE static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) { // FIXME make sure lhs is sequentially stored const int cols = dest.cols(); for (int j=0; j struct ei_outer_product_selector { template<> struct ei_outer_product_selector { template - static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) { + EIGEN_DONT_INLINE static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) { // FIXME make sure rhs is sequentially stored const int rows = dest.rows(); for (int i=0; isetZero(); } +template +void MatrixBase::makeHouseholderInPlace(RealScalar *tau, Scalar *beta) +{ + VectorBlock::ret> essentialPart(derived(), 1, size()-1); + makeHouseholder(&essentialPart, tau, beta); +} + +/** Computes the elementary reflector H such that: + * \f$ H *this = [ beta 0 ... 0]^T \f$ + * where the transformation H is: + * \f$ H = I - tau v v^*\f$ + * and the vector v is: + * \f$ v^T = [1 essential^T] \f$ + * + * On output: + * \param essential the essential part of the vector \c v + * \param tau the scaling factor of the householder transformation + * \param beta the result of H * \c *this + * + * \sa MatrixBase::makeHouseholderInPlace(), MatrixBase::applyHouseholderOnTheLeft(), + * MatrixBase::applyHouseholderOnTheRight() + */ template template void MatrixBase::makeHouseholder( EssentialPart *essential, - RealScalar *beta) const + RealScalar *tau, + Scalar *beta) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) RealScalar _squaredNorm = squaredNorm(); @@ -62,35 +85,38 @@ void MatrixBase::makeHouseholder( VectorBlock tail(derived(), 1, size()-1); *essential = tail / c0; const RealScalar c0abs2 = ei_abs2(c0); - *beta = RealScalar(2) * c0abs2 / (c0abs2 + _squaredNorm - ei_abs2(coeff(0))); + *tau = RealScalar(2) * c0abs2 / (c0abs2 + _squaredNorm - ei_abs2(coeff(0))); + *beta = coeff(0) - c0; } template template void MatrixBase::applyHouseholderOnTheLeft( const EssentialPart& essential, - const RealScalar& beta) + const RealScalar& tau, + Scalar* workspace) { - Matrix tmp(cols()); + Map > tmp(workspace,cols()); Block bottom(derived(), 1, 0, rows()-1, cols()); - tmp = (essential.adjoint() * bottom).lazy(); + tmp.noalias() = essential.adjoint() * bottom; tmp += row(0); - row(0) -= beta * tmp; - bottom -= beta * essential * tmp; + row(0) -= tau * tmp; + bottom.noalias() -= tau * essential * tmp; } template template void MatrixBase::applyHouseholderOnTheRight( const EssentialPart& essential, - const RealScalar& beta) + const RealScalar& tau, + Scalar* workspace) { - Matrix tmp(rows()); + Map > tmp(workspace,rows()); Block right(derived(), 0, 1, rows(), cols()-1); - tmp = (right * essential.conjugate()).lazy(); + tmp.noalias() = right * essential.conjugate(); tmp += col(0); - col(0) -= beta * tmp; - right -= beta * tmp * essential.transpose(); + col(0) -= tau * tmp; + right.noalias() -= tau * tmp * essential.transpose(); } #endif // EIGEN_HOUSEHOLDER_H diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index d37a019ff..0c6142129 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -45,11 +45,15 @@ template class HouseholderQR { public: + enum { + MinSizeAtCompileTime = EIGEN_ENUM_MIN(MatrixType::ColsAtCompileTime,MatrixType::RowsAtCompileTime) + }; + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Block MatrixRBlockType; typedef Matrix MatrixTypeR; - typedef Matrix VectorType; + typedef Matrix HCoeffsType; /** * \brief Default Constructor. @@ -61,7 +65,7 @@ template class HouseholderQR HouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), - m_hCoeffs(matrix.cols()), + m_hCoeffs(std::min(matrix.rows(),matrix.cols())), m_isInitialized(false) { compute(matrix); @@ -105,7 +109,7 @@ template class HouseholderQR protected: MatrixType m_qr; - VectorType m_hCoeffs; + HCoeffsType m_hCoeffs; bool m_isInitialized; }; @@ -114,65 +118,25 @@ template class HouseholderQR template HouseholderQR& HouseholderQR::compute(const MatrixType& matrix) { - m_qr = matrix; - m_hCoeffs.resize(matrix.cols()); - int rows = matrix.rows(); int cols = matrix.cols(); - RealScalar eps2 = precision()*precision(); + int size = std::min(rows,cols); + + m_qr = matrix; + m_hCoeffs.resize(size); + Matrix temp(cols); - for (int k = 0; k < cols; ++k) + for (int k = 0; k < size; ++k) { - int remainingSize = rows-k; + int remainingRows = rows - k; + int remainingCols = cols - k -1; - RealScalar beta; - Scalar v0 = m_qr.col(k).coeff(k); - - if (remainingSize==1) - { - if (NumTraits::IsComplex) - { - // Householder transformation on the remaining single scalar - beta = ei_abs(v0); - if (ei_real(v0)>0) - beta = -beta; - m_qr.coeffRef(k,k) = beta; - m_hCoeffs.coeffRef(k) = (beta - v0) / beta; - } - else - { - m_hCoeffs.coeffRef(k) = 0; - } - } - else if ((beta=m_qr.col(k).end(remainingSize-1).squaredNorm())>eps2) - // FIXME what about ei_imag(v0) ?? - { - // form k-th Householder vector - beta = ei_sqrt(ei_abs2(v0)+beta); - if (ei_real(v0)>=0.) - beta = -beta; - m_qr.col(k).end(remainingSize-1) /= v0-beta; - m_qr.coeffRef(k,k) = beta; - Scalar h = m_hCoeffs.coeffRef(k) = (beta - v0) / beta; - - // apply the Householder transformation (I - h v v') to remaining columns, i.e., - // R <- (I - h v v') * R where v = [1,m_qr(k+1,k), m_qr(k+2,k), ...] - int remainingCols = cols - k -1; - if (remainingCols>0) - { - m_qr.coeffRef(k,k) = Scalar(1); - temp.end(remainingCols) = (m_qr.col(k).end(remainingSize).adjoint() - * m_qr.corner(BottomRight, remainingSize, remainingCols)).lazy(); - m_qr.corner(BottomRight, remainingSize, remainingCols) -= (ei_conj(h) * m_qr.col(k).end(remainingSize) - * temp.end(remainingCols)).lazy(); - m_qr.coeffRef(k,k) = beta; - } - } - else - { - m_hCoeffs.coeffRef(k) = 0; - } + m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &m_qr.coeffRef(k,k)); + + if (remainingRows>1 && remainingCols>0) + m_qr.corner(BottomRight, remainingRows, remainingCols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); } m_isInitialized = true; return *this; @@ -187,12 +151,20 @@ void HouseholderQR::solve( { ei_assert(m_isInitialized && "HouseholderQR is not initialized."); const int rows = m_qr.rows(); + const int cols = b.cols(); ei_assert(b.rows() == rows); - result->resize(rows, b.cols()); + result->resize(rows, cols); - // TODO(keir): There is almost certainly a faster way to multiply by - // Q^T without explicitly forming matrixQ(). Investigate. - *result = matrixQ().transpose()*b; + *result = b; + + Matrix temp(cols); + for (int k = 0; k < cols; ++k) + { + int remainingSize = rows-k; + + result->corner(BottomRight, remainingSize, cols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(0)); + } const int rank = std::min(result->rows(), result->cols()); m_qr.corner(TopLeft, rank, rank) @@ -214,15 +186,10 @@ MatrixType HouseholderQR::matrixQ() const Matrix temp(cols); for (int k = cols-1; k >= 0; k--) { - // to make easier the computation of the transformation, let's temporarily - // overwrite m_qr(k,k) such that the end of m_qr.col(k) is exactly our Householder vector. - Scalar beta = m_qr.coeff(k,k); - m_qr.const_cast_derived().coeffRef(k,k) = 1; - int endLength = rows-k; + int remainingSize = rows-k; - temp.end(cols-k) = (m_qr.col(k).end(endLength).adjoint() * res.corner(BottomRight,endLength, cols-k)).lazy(); - res.corner(BottomRight,endLength, cols-k) -= ((m_hCoeffs.coeff(k) * m_qr.col(k).end(endLength)) * temp.end(cols-k)).lazy(); - m_qr.const_cast_derived().coeffRef(k,k) = beta; + res.corner(BottomRight, remainingSize, cols-k) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); } return res; } diff --git a/test/qr.cpp b/test/qr.cpp index 88a447c4b..99df1d89b 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -37,8 +37,8 @@ template void qr(const MatrixType& m) MatrixType a = MatrixType::Random(rows,cols); HouseholderQR qrOfA(a); - VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR()); - VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR()); + VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR().toDense()); + VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR().toDense()); SquareMatrixType b = a.adjoint() * a; @@ -59,7 +59,7 @@ template void qr_invertible() { /* this test covers the following files: QR.h */ typedef typename NumTraits::Real RealScalar; - int size = ei_random(10,200); + int size = ei_random(10,50); MatrixType m1(size, size), m2(size, size), m3(size, size); m1 = MatrixType::Random(size,size); @@ -74,7 +74,6 @@ template void qr_invertible() HouseholderQR qr(m1); m3 = MatrixType::Random(size,size); qr.solve(m3, &m2); - //std::cerr << m3 - m1*m2 << "\n\n"; VERIFY_IS_APPROX(m3, m1*m2); } @@ -91,20 +90,18 @@ template void qr_verify_assert() void test_qr() { for(int i = 0; i < 1; i++) { + // FIXME : very weird bug here // CALL_SUBTEST( qr(Matrix2f()) ); -// CALL_SUBTEST( qr(Matrix4d()) ); -// CALL_SUBTEST( qr(MatrixXf(12,8)) ); -// CALL_SUBTEST( qr(MatrixXcd(5,5)) ); -// CALL_SUBTEST( qr(MatrixXcd(7,3)) ); - CALL_SUBTEST( qr(MatrixXf(47,47)) ); + CALL_SUBTEST( qr(Matrix4d()) ); + CALL_SUBTEST( qr(MatrixXcd(17,7)) ); + CALL_SUBTEST( qr(MatrixXf(47,40)) ); } for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST( qr_invertible() ); CALL_SUBTEST( qr_invertible() ); - // TODO fix issue with complex -// CALL_SUBTEST( qr_invertible() ); -// CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); } CALL_SUBTEST(qr_verify_assert()); From e125c199bbe3c0b61c8732c7603b66745c4582fe Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 17 Aug 2009 09:16:41 +0200 Subject: [PATCH 54/58] add EIGEN_TRANSFORM_PLUGIN --- Eigen/src/Geometry/Transform.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 642f307ed..9eb8ed535 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -460,6 +460,11 @@ public: inline const Block translationExt() const { return m_matrix.template block(0,Dim); } + + #ifdef EIGEN_TRANSFORM_PLUGIN + #include EIGEN_TRANSFORM_PLUGIN + #endif + }; /** \ingroup Geometry_Module */ From ff0f005d4c0bd46e88d050b9f147eab810f4814d Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 17 Aug 2009 17:04:32 +0200 Subject: [PATCH 55/58] change the make householder algorithm so that the remaining coefficient is real, and make Tridiagonalization use it --- Eigen/src/Core/MatrixBase.h | 8 +-- Eigen/src/Householder/Householder.h | 34 ++++++------ Eigen/src/QR/QR.h | 23 +++++---- Eigen/src/QR/Tridiagonalization.h | 80 +++++++---------------------- test/eigensolver_selfadjoint.cpp | 2 +- test/householder.cpp | 31 +++++++---- test/qr.cpp | 2 +- 7 files changed, 75 insertions(+), 105 deletions(-) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 9e92c043f..1f4c6bf7a 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -786,17 +786,17 @@ template class MatrixBase ////////// Householder module /////////// - void makeHouseholderInPlace(RealScalar *tau, Scalar *beta); + void makeHouseholderInPlace(Scalar *tau, RealScalar *beta); template void makeHouseholder(EssentialPart *essential, - RealScalar *tau, Scalar *beta) const; + Scalar *tau, RealScalar *beta) const; template void applyHouseholderOnTheLeft(const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace); template void applyHouseholderOnTheRight(const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace); ///////// Jacobi module ///////// diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index a7bbe96ce..769ba3d43 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -42,7 +42,7 @@ void makeTrivialHouseholder( } template -void MatrixBase::makeHouseholderInPlace(RealScalar *tau, Scalar *beta) +void MatrixBase::makeHouseholderInPlace(Scalar *tau, RealScalar *beta) { VectorBlock::ret> essentialPart(derived(), 1, size()-1); makeHouseholder(&essentialPart, tau, beta); @@ -67,33 +67,35 @@ template template void MatrixBase::makeHouseholder( EssentialPart *essential, - RealScalar *tau, - Scalar *beta) const + Scalar *tau, + RealScalar *beta) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) - RealScalar _squaredNorm = squaredNorm(); - Scalar c0; - if(ei_abs2(coeff(0)) <= ei_abs2(precision()) * _squaredNorm) + VectorBlock tail(derived(), 1, size()-1); + + RealScalar tailSqNorm; + Scalar c0 = coeff(0); + + if( (size()==1 || (tailSqNorm=tail.squaredNorm()) == RealScalar(0)) && ei_imag(c0)==RealScalar(0)) { - c0 = ei_sqrt(_squaredNorm); + *tau = 0; + *beta = ei_real(c0); } else { - Scalar sign = coeff(0) / ei_abs(coeff(0)); - c0 = coeff(0) + sign * ei_sqrt(_squaredNorm); + *beta = ei_sqrt(ei_abs2(c0) + tailSqNorm); + if (ei_real(c0)>=0.) + *beta = -*beta; + *essential = tail / (c0 - *beta); + *tau = ei_conj((*beta - c0) / *beta); } - VectorBlock tail(derived(), 1, size()-1); - *essential = tail / c0; - const RealScalar c0abs2 = ei_abs2(c0); - *tau = RealScalar(2) * c0abs2 / (c0abs2 + _squaredNorm - ei_abs2(coeff(0))); - *beta = coeff(0) - c0; } template template void MatrixBase::applyHouseholderOnTheLeft( const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace) { Map > tmp(workspace,cols()); @@ -108,7 +110,7 @@ template template void MatrixBase::applyHouseholderOnTheRight( const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace) { Map > tmp(workspace,rows()); diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 0c6142129..90e6f8132 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -53,7 +53,7 @@ template class HouseholderQR typedef typename MatrixType::RealScalar RealScalar; typedef Block MatrixRBlockType; typedef Matrix MatrixTypeR; - typedef Matrix HCoeffsType; + typedef Matrix HCoeffsType; /** * \brief Default Constructor. @@ -132,11 +132,13 @@ HouseholderQR& HouseholderQR::compute(const MatrixType& int remainingRows = rows - k; int remainingCols = cols - k -1; - m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &m_qr.coeffRef(k,k)); - - if (remainingRows>1 && remainingCols>0) - m_qr.corner(BottomRight, remainingRows, remainingCols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + RealScalar beta; + m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); + m_qr.coeffRef(k,k) = beta; + + // apply H to remaining part of m_qr from the left + m_qr.corner(BottomRight, remainingRows, remainingCols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); } m_isInitialized = true; return *this; @@ -163,7 +165,7 @@ void HouseholderQR::solve( int remainingSize = rows-k; result->corner(BottomRight, remainingSize, cols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(0)); + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); } const int rank = std::min(result->rows(), result->cols()); @@ -177,8 +179,8 @@ template MatrixType HouseholderQR::matrixQ() const { ei_assert(m_isInitialized && "HouseholderQR is not initialized."); - // compute the product Q_0 Q_1 ... Q_n-1, - // where Q_k is the k-th Householder transformation I - h_k v_k v_k' + // compute the product H'_0 H'_1 ... H'_n-1, + // where H_k is the k-th Householder transformation I - h_k v_k v_k' // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] int rows = m_qr.rows(); int cols = m_qr.cols(); @@ -187,9 +189,8 @@ MatrixType HouseholderQR::matrixQ() const for (int k = cols-1; k >= 0; k--) { int remainingSize = rows-k; - res.corner(BottomRight, remainingSize, cols-k) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); } return res; } diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index d98322fac..4cc7433c1 100644 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -198,65 +198,29 @@ void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& { assert(matA.rows()==matA.cols()); int n = matA.rows(); - for (int i = 0; i aux(n); + for (int i = 0; i(1))) - { - hCoeffs.coeffRef(i) = 0.; - } - else - { - Scalar v0 = matA.col(i).coeff(i+1); - RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2); - if (ei_real(v0)>=0.) - beta = -beta; - matA.col(i).end(n-(i+2)) *= (Scalar(1)/(v0-beta)); - matA.col(i).coeffRef(i+1) = beta; - Scalar h = (beta - v0) / beta; - // end of the householder transformation + hCoeffs.end(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView() + * (ei_conj(h) * matA.col(i).end(remainingSize))); - // Apply similarity transformation to remaining columns, - // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1) - matA.col(i).coeffRef(i+1) = 1; + hCoeffs.end(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.end(remainingSize).dot(matA.col(i).end(remainingSize)))) * matA.col(i).end(n-i-1); - hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template selfadjointView() - * (h * matA.col(i).end(n-i-1))); + matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView() + .rankUpdate(matA.col(i).end(remainingSize), hCoeffs.end(remainingSize), -1); - hCoeffs.end(n-i-1) += (h*Scalar(-0.5)*(hCoeffs.end(n-i-1).dot(matA.col(i).end(n-i-1)))) * matA.col(i).end(n-i-1); - - matA.corner(BottomRight, n-i-1, n-i-1).template selfadjointView() - .rankUpdate(matA.col(i).end(n-i-1), hCoeffs.end(n-i-1), -1); - - // note: at that point matA(i+1,i+1) is the (i+1)-th element of the final diagonal - // note: the sequence of the beta values leads to the subdiagonal entries - matA.col(i).coeffRef(i+1) = beta; - - hCoeffs.coeffRef(i) = h; - } - } - if (NumTraits::IsComplex) - { - // Householder transformation on the remaining single scalar - int i = n-2; - Scalar v0 = matA.col(i).coeff(i+1); - RealScalar beta = ei_abs(v0); - if (ei_real(v0)>=RealScalar(0)) - beta = -beta; matA.col(i).coeffRef(i+1) = beta; - // FIXME comparing against 1 - if(ei_isMuchSmallerThan(beta, Scalar(1))) hCoeffs.coeffRef(i) = Scalar(0); - else hCoeffs.coeffRef(i) = (beta - v0) / beta; - } - else - { - hCoeffs.coeffRef(n-2) = 0; + hCoeffs.coeffRef(i) = h; } } @@ -280,16 +244,8 @@ void Tridiagonalization::matrixQInPlace(MatrixBase* q) con Matrix aux(n); for (int i = n-2; i>=0; i--) { - Scalar tmp = m_matrix.coeff(i+1,i); - m_matrix.const_cast_derived().coeffRef(i+1,i) = 1; - - aux.end(n-i-1) = (m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy(); - // rank one update, TODO ! make it works efficiently as expected - for (int j=i+1;j void householder(const MatrixType& m) { + static bool even = true; + even = !even; /* this test covers the following files: Householder.h */ @@ -38,46 +40,55 @@ template void householder(const MatrixType& m) typedef Matrix VectorType; typedef Matrix::ret, 1> EssentialVectorType; typedef Matrix SquareMatrixType; + + Matrix _tmp(std::max(rows,cols)); + Scalar* tmp = &_tmp.coeffRef(0,0); - RealScalar beta; + Scalar beta; + RealScalar alpha; EssentialVectorType essential; VectorType v1 = VectorType::Random(rows), v2; v2 = v1; - v1.makeHouseholder(&essential, &beta); - v1.applyHouseholderOnTheLeft(essential,beta); - + v1.makeHouseholder(&essential, &beta, &alpha); + v1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(v1.norm(), v2.norm()); VERIFY_IS_MUCH_SMALLER_THAN(v1.end(rows-1).norm(), v1.norm()); v1 = VectorType::Random(rows); v2 = v1; - v1.applyHouseholderOnTheLeft(essential,beta); + v1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(v1.norm(), v2.norm()); MatrixType m1(rows, cols), m2(rows, cols); v1 = VectorType::Random(rows); + if(even) v1.end(rows-1).setZero(); m1.colwise() = v1; m2 = m1; - m1.col(0).makeHouseholder(&essential, &beta); - m1.applyHouseholderOnTheLeft(essential,beta); + m1.col(0).makeHouseholder(&essential, &beta, &alpha); + m1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(m1.norm(), m2.norm()); VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm()); + VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m1(0,0)), ei_real(m1(0,0))); + VERIFY_IS_APPROX(ei_real(m1(0,0)), alpha); v1 = VectorType::Random(rows); + if(even) v1.end(rows-1).setZero(); SquareMatrixType m3(rows,rows), m4(rows,rows); m3.rowwise() = v1.transpose(); m4 = m3; - m3.row(0).makeHouseholder(&essential, &beta); - m3.applyHouseholderOnTheRight(essential,beta); + m3.row(0).makeHouseholder(&essential, &beta, &alpha); + m3.applyHouseholderOnTheRight(essential,beta,tmp); VERIFY_IS_APPROX(m3.norm(), m4.norm()); VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm()); + VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0))); + VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha); } void test_householder() { - for(int i = 0; i < g_repeat; i++) { + for(int i = 0; i < 2*g_repeat; i++) { CALL_SUBTEST( householder(Matrix()) ); CALL_SUBTEST( householder(Matrix()) ); CALL_SUBTEST( householder(Matrix()) ); diff --git a/test/qr.cpp b/test/qr.cpp index 99df1d89b..f004a36ca 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -93,8 +93,8 @@ void test_qr() // FIXME : very weird bug here // CALL_SUBTEST( qr(Matrix2f()) ); CALL_SUBTEST( qr(Matrix4d()) ); - CALL_SUBTEST( qr(MatrixXcd(17,7)) ); CALL_SUBTEST( qr(MatrixXf(47,40)) ); + CALL_SUBTEST( qr(MatrixXcd(17,7)) ); } for(int i = 0; i < g_repeat; i++) { From d56be9c128dbc88821d2ef4fa7d48ced72923d17 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 17 Aug 2009 17:41:01 +0200 Subject: [PATCH 56/58] * make HessenbergDecomposition uses the Householder module * bugfix in ei_blas_traits for .conjugate().conjugate() --- Eigen/src/Core/util/BlasUtil.h | 2 +- Eigen/src/QR/HessenbergDecomposition.h | 91 +++++++------------------- Eigen/src/QR/Tridiagonalization.h | 1 - 3 files changed, 24 insertions(+), 70 deletions(-) diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h index f4690c0ca..94154108c 100644 --- a/Eigen/src/Core/util/BlasUtil.h +++ b/Eigen/src/Core/util/BlasUtil.h @@ -182,7 +182,7 @@ struct ei_blas_traits, NestedXpr> > enum { IsComplex = NumTraits::IsComplex, - NeedToConjugate = IsComplex + NeedToConjugate = Base::NeedToConjugate ? 0 : IsComplex }; static inline ExtractType extract(const XprType& x) { return Base::extract(x._expression()); } static inline Scalar extractScalarFactor(const XprType& x) { return ei_conj(Base::extractScalarFactor(x._expression())); } diff --git a/Eigen/src/QR/HessenbergDecomposition.h b/Eigen/src/QR/HessenbergDecomposition.h index 6ba90fb3a..6b261a8a7 100644 --- a/Eigen/src/QR/HessenbergDecomposition.h +++ b/Eigen/src/QR/HessenbergDecomposition.h @@ -93,7 +93,7 @@ template class HessenbergDecomposition * * \sa packedMatrix() */ - CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; } + CoeffVectorType householderCoefficients() const { return m_hCoeffs; } /** \returns the internal result of the decomposition. * @@ -112,8 +112,8 @@ template class HessenbergDecomposition */ const MatrixType& packedMatrix(void) const { return m_matrix; } - MatrixType matrixQ(void) const; - MatrixType matrixH(void) const; + MatrixType matrixQ() const; + MatrixType matrixH() const; private: @@ -143,87 +143,42 @@ void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVector { assert(matA.rows()==matA.cols()); int n = matA.rows(); - for (int i = 0; i temp(n); + for (int i = 0; i(1))) - { - hCoeffs.coeffRef(i) = 0.; - } - else - { - Scalar v0 = matA.col(i).coeff(i+1); - RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2); - if (ei_real(v0)>=0.) - beta = -beta; - matA.col(i).end(n-(i+2)) *= (Scalar(1)/(v0-beta)); - matA.col(i).coeffRef(i+1) = beta; - Scalar h = (beta - v0) / beta; - // end of the householder transformation - - // Apply similarity transformation to remaining columns, - // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1) - matA.col(i).coeffRef(i+1) = 1; - - // first let's do A = H' A - matA.corner(BottomRight,n-i-1,n-i-1) -= ((ei_conj(h) * matA.col(i).end(n-i-1)) * - (matA.col(i).end(n-i-1).adjoint() * matA.corner(BottomRight,n-i-1,n-i-1))).lazy(); - - // now let's do A = A H - matA.corner(BottomRight,n,n-i-1) -= ((matA.corner(BottomRight,n,n-i-1) * matA.col(i).end(n-i-1)) - * (h * matA.col(i).end(n-i-1).adjoint())).lazy(); - - matA.col(i).coeffRef(i+1) = beta; - hCoeffs.coeffRef(i) = h; - } - } - if (NumTraits::IsComplex) - { - // Householder transformation on the remaining single scalar - int i = n-2; - Scalar v0 = matA.coeff(i+1,i); - - RealScalar beta = ei_sqrt(ei_abs2(v0)); - if (ei_real(v0)>=0.) - beta = -beta; - Scalar h = (beta - v0) / beta; + int remainingSize = n-i-1; + RealScalar beta; + Scalar h; + matA.col(i).end(remainingSize).makeHouseholderInPlace(&h, &beta); + matA.col(i).coeffRef(i+1) = beta; hCoeffs.coeffRef(i) = h; - // A = H* A - matA.corner(BottomRight,n-i-1,n-i) -= ei_conj(h) * matA.corner(BottomRight,n-i-1,n-i); + // Apply similarity transformation to remaining columns, + // i.e., compute A = H A H' + + // A = H A + matA.corner(BottomRight, remainingSize, remainingSize) + .applyHouseholderOnTheLeft(matA.col(i).end(remainingSize-1), h, &temp.coeffRef(0)); - // A = A H - matA.col(n-1) -= h * matA.col(n-1); - } - else - { - hCoeffs.coeffRef(n-2) = 0; + // A = A H' + matA.corner(BottomRight, n, remainingSize) + .applyHouseholderOnTheRight(matA.col(i).end(remainingSize-1).conjugate(), ei_conj(h), &temp.coeffRef(0)); } } /** reconstructs and returns the matrix Q */ template typename HessenbergDecomposition::MatrixType -HessenbergDecomposition::matrixQ(void) const +HessenbergDecomposition::matrixQ() const { int n = m_matrix.rows(); MatrixType matQ = MatrixType::Identity(n,n); Matrix temp(n); for (int i = n-2; i>=0; i--) { - Scalar tmp = m_matrix.coeff(i+1,i); - m_matrix.const_cast_derived().coeffRef(i+1,i) = 1; - - int rs = n-i-1; - temp.end(rs) = (m_matrix.col(i).end(rs).adjoint() * matQ.corner(BottomRight,rs,rs)).lazy(); - matQ.corner(BottomRight,rs,rs) -= (m_hCoeffs.coeff(i) * m_matrix.col(i).end(rs) * temp.end(rs)).lazy(); - - m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp; + matQ.corner(BottomRight,n-i-1,n-i-1) + .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0)); } return matQ; } @@ -237,7 +192,7 @@ HessenbergDecomposition::matrixQ(void) const */ template typename HessenbergDecomposition::MatrixType -HessenbergDecomposition::matrixH(void) const +HessenbergDecomposition::matrixH() const { // FIXME should this function (and other similar) rather take a matrix as argument // and fill it (to avoid temporaries) diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index 4cc7433c1..2053505f6 100644 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -201,7 +201,6 @@ void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& Matrix aux(n); for (int i = 0; i Date: Tue, 18 Aug 2009 11:09:20 +0100 Subject: [PATCH 57/58] Correct syntax error in doxygen comment. --- unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h index cbe35a110..c2b77636c 100644 --- a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h +++ b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h @@ -108,7 +108,7 @@ void ei_pseudo_inverse(const CMatrix &C, CINVMatrix &CINV) /** \ingroup IterativeSolvers_Module * Constrained conjugate gradient * - * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the contraint \f$ Cx <= f @\$ + * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the contraint \f$ Cx \le f \f$ */ template From 59a0c4a0d2727016814bf62a2b8fbc8b8420fd32 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Tue, 18 Aug 2009 15:30:38 +0100 Subject: [PATCH 58/58] Add new unsupported modules to doc/unsupported_modules.dox --- doc/unsupported_modules.dox | 16 ++++++++++++++-- unsupported/Eigen/MoreVectorization | 2 +- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/doc/unsupported_modules.dox b/doc/unsupported_modules.dox index 61ea9ada2..856b05831 100644 --- a/doc/unsupported_modules.dox +++ b/doc/unsupported_modules.dox @@ -14,9 +14,21 @@ namespace Eigen { * \defgroup AdolcForward_Module Adolc forward module */ /** \ingroup Unsupported_modules - * \defgroup IterativeSolvers_Module Iterative solvers module */ + * \defgroup AlignedVector3_Module Aligned vector3 module */ + +/** \ingroup Unsupported_modules + * \defgroup AutoDiff_Module Auto Diff module */ /** \ingroup Unsupported_modules * \defgroup BVH_Module BVH module */ -} // namespace Eigen \ No newline at end of file +/** \ingroup Unsupported_modules + * \defgroup IterativeSolvers_Module Iterative solvers module */ + +/** \ingroup Unsupported_modules + * \defgroup MatrixFunctions_Module Matrix functions module */ + +/** \ingroup Unsupported_modules + * \defgroup MoreVectorization More vectorization module */ + +} // namespace Eigen diff --git a/unsupported/Eigen/MoreVectorization b/unsupported/Eigen/MoreVectorization index e60526a3d..26a01cd29 100644 --- a/unsupported/Eigen/MoreVectorization +++ b/unsupported/Eigen/MoreVectorization @@ -6,7 +6,7 @@ namespace Eigen { /** \ingroup Unsupported_modules - * \defgroup MoreVectorization additional vectorization module + * \defgroup MoreVectorization More vectorization module */ #include "src/MoreVectorization/MathFunctions.h"