From 5960befc206ac7405841a9da2436f377e4df694f Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 19 Feb 2014 21:42:29 +0100 Subject: [PATCH 01/24] More int versus Index fixes --- Eigen/src/SparseCore/SparseMatrix.h | 9 +++++---- Eigen/src/SparseCore/SparsePermutation.h | 4 ++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 0e7e760fb..9ac18bcf6 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -223,7 +223,7 @@ class SparseMatrix if(isCompressed()) { - reserve(VectorXi::Constant(outerSize(), 2)); + reserve(Matrix::Constant(outerSize(), 2)); } return insertUncompressed(row,col); } @@ -939,12 +939,13 @@ void set_from_triplets(const InputIterator& begin, const InputIterator& end, Spa EIGEN_UNUSED_VARIABLE(Options); enum { IsRowMajor = SparseMatrixType::IsRowMajor }; typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::Index Index; SparseMatrix trMat(mat.rows(),mat.cols()); if(begin!=end) { // pass 1: count the nnz per inner-vector - VectorXi wi(trMat.outerSize()); + Matrix wi(trMat.outerSize()); wi.setZero(); for(InputIterator it(begin); it!=end; ++it) { @@ -1018,7 +1019,7 @@ void SparseMatrix::sumupDuplicates() { eigen_assert(!isCompressed()); // TODO, in practice we should be able to use m_innerNonZeros for that task - VectorXi wi(innerSize()); + Matrix wi(innerSize()); wi.fill(-1); Index count = 0; // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers @@ -1081,7 +1082,7 @@ EIGEN_DONT_INLINE SparseMatrix& SparseMatrix positions(dest.outerSize()); for (Index j=0; j tmp(m_matrix.rows(), m_matrix.cols()); - VectorXi sizes(m_matrix.outerSize()); + Matrix sizes(m_matrix.outerSize()); for(Index j=0; j tmp(m_matrix.rows(), m_matrix.cols()); - VectorXi sizes(tmp.outerSize()); + Matrix sizes(tmp.outerSize()); sizes.setZero(); PermutationMatrix perm; if((Side==OnTheLeft) ^ Transposed) From 3e439889e0dd10f67d328f3c688178d1e6c091d8 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Mon, 24 Feb 2014 13:12:10 +0100 Subject: [PATCH 02/24] Specify what non-resizeable objects are in transposeInPlace and adjointInPlace (cf bug #749) --- Eigen/src/Core/Transpose.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 976708a0f..aba3f6670 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -293,7 +293,8 @@ struct inplace_transpose_selector { // non square matrix * Notice however that this method is only useful if you want to replace a matrix by its own transpose. * If you just need the transpose of a matrix, use transpose(). * - * \note if the matrix is not square, then \c *this must be a resizable matrix. + * \note if the matrix is not square, then \c *this must be a resizable matrix. + * This excludes (non-square) fixed-size matrices, block-expressions and maps. * * \sa transpose(), adjoint(), adjointInPlace() */ template @@ -324,6 +325,7 @@ inline void DenseBase::transposeInPlace() * If you just need the adjoint of a matrix, use adjoint(). * * \note if the matrix is not square, then \c *this must be a resizable matrix. + * This excludes (non-square) fixed-size matrices, block-expressions and maps. * * \sa transpose(), adjoint(), transposeInPlace() */ template From 6fecb6f1b6d17883525f3f34dd5c1860faeecbd8 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Mon, 24 Feb 2014 14:10:17 +0000 Subject: [PATCH 03/24] Fix bug #748 - array_5 test fails for seed 1392781168. --- test/array.cpp | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/test/array.cpp b/test/array.cpp index f1deda7e3..5f49fc1ea 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -173,21 +173,14 @@ template void array_real(const ArrayType& m) Scalar s1 = internal::random(); // these tests are mostly to check possible compilation issues. -// VERIFY_IS_APPROX(m1.sin(), std::sin(m1)); VERIFY_IS_APPROX(m1.sin(), sin(m1)); -// VERIFY_IS_APPROX(m1.cos(), std::cos(m1)); VERIFY_IS_APPROX(m1.cos(), cos(m1)); -// VERIFY_IS_APPROX(m1.asin(), std::asin(m1)); VERIFY_IS_APPROX(m1.asin(), asin(m1)); -// VERIFY_IS_APPROX(m1.acos(), std::acos(m1)); VERIFY_IS_APPROX(m1.acos(), acos(m1)); -// VERIFY_IS_APPROX(m1.tan(), std::tan(m1)); VERIFY_IS_APPROX(m1.tan(), tan(m1)); VERIFY_IS_APPROX(cos(m1+RealScalar(3)*m2), cos((m1+RealScalar(3)*m2).eval())); -// VERIFY_IS_APPROX(std::cos(m1+RealScalar(3)*m2), std::cos((m1+RealScalar(3)*m2).eval())); -// VERIFY_IS_APPROX(m1.abs().sqrt(), std::sqrt(std::abs(m1))); VERIFY_IS_APPROX(m1.abs().sqrt(), sqrt(abs(m1))); VERIFY_IS_APPROX(m1.abs(), sqrt(numext::abs2(m1))); @@ -196,9 +189,10 @@ template void array_real(const ArrayType& m) if(!NumTraits::IsComplex) VERIFY_IS_APPROX(numext::real(m1), m1); - VERIFY_IS_APPROX(m1.abs().log() , log(abs(m1))); + // shift argument of logarithm so that it is not zero + Scalar smallNumber = NumTraits::dummy_precision(); + VERIFY_IS_APPROX((m1.abs() + smallNumber).log() , log(abs(m1) + smallNumber)); -// VERIFY_IS_APPROX(m1.exp(), std::exp(m1)); VERIFY_IS_APPROX(m1.exp() * m2.exp(), exp(m1+m2)); VERIFY_IS_APPROX(m1.exp(), exp(m1)); VERIFY_IS_APPROX(m1.exp() / m2.exp(),(m1-m2).exp()); @@ -242,7 +236,6 @@ template void array_complex(const ArrayType& m) m2(i,j) = sqrt(m1(i,j)); VERIFY_IS_APPROX(m1.sqrt(), m2); -// VERIFY_IS_APPROX(m1.sqrt(), std::sqrt(m1)); VERIFY_IS_APPROX(m1.sqrt(), Eigen::sqrt(m1)); } From 21fecd5252061c05cfb9f91e4e9becf25089949d Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 24 Feb 2014 17:12:17 +0100 Subject: [PATCH 04/24] Workaround clang ABI change with unsed arguments (ugly fix) --- .../Core/products/GeneralBlockPanelKernel.h | 12 +++++++ Eigen/src/Core/products/GeneralMatrixVector.h | 31 ++++++++++++------- 2 files changed, 31 insertions(+), 12 deletions(-) diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 780fa74d3..10cd34e9e 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -1130,6 +1130,10 @@ EIGEN_DONT_INLINE void gemm_pack_lhs=depth && offset<=stride)); eigen_assert( (StorageOrder==RowMajor) || ((Pack1%PacketSize)==0 && Pack1<=4*PacketSize) ); +#ifdef __clang__ + // Workaround clang ABI change with unsed arguments + if(!PanelMode) depth += stride + offset; +#endif conj_if::IsComplex && Conjugate> cj; const_blas_data_mapper lhs(_lhs,lhsStride); Index count = 0; @@ -1216,6 +1220,10 @@ EIGEN_DONT_INLINE void gemm_pack_rhs=depth && offset<=stride)); +#ifdef __clang__ + // Workaround clang ABI change with unsed arguments + if(!PanelMode) depth += stride + offset; +#endif conj_if::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; Index count = 0; @@ -1267,6 +1275,10 @@ EIGEN_DONT_INLINE void gemm_pack_rhs=depth && offset<=stride)); +#ifdef __clang__ + // Workaround clang ABI change with unsed arguments + if(!PanelMode) depth += stride + offset; +#endif conj_if::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; Index count = 0; diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index c8697c913..c3e1e7142 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h @@ -80,11 +80,8 @@ EIGEN_DONT_INLINE static void run( Index rows, Index cols, const LhsScalar* lhs, Index lhsStride, const RhsScalar* rhs, Index rhsIncr, - ResScalar* res, Index - #ifdef EIGEN_INTERNAL_DEBUGGING - resIncr - #endif - , RhsScalar alpha); + ResScalar* res, Index resIncr, + RhsScalar alpha); }; template @@ -92,13 +89,17 @@ EIGEN_DONT_INLINE void general_matrix_vector_product Date: Mon, 24 Feb 2014 18:13:49 +0100 Subject: [PATCH 05/24] Implement bug #317: use a template function call to suppress unused variable warnings. This also fix the issue of the previous changeset in a much nicer way. --- .../Core/products/GeneralBlockPanelKernel.h | 18 ++++++------------ Eigen/src/Core/products/GeneralMatrixVector.h | 15 ++------------- Eigen/src/Core/util/Macros.h | 7 ++++++- 3 files changed, 14 insertions(+), 26 deletions(-) diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 10cd34e9e..08cc14bd7 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -1128,12 +1128,10 @@ EIGEN_DONT_INLINE void gemm_pack_lhs::size }; EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK LHS"); + EIGEN_UNUSED_VARIABLE(stride); + EIGEN_UNUSED_VARIABLE(offset); eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride)); eigen_assert( (StorageOrder==RowMajor) || ((Pack1%PacketSize)==0 && Pack1<=4*PacketSize) ); -#ifdef __clang__ - // Workaround clang ABI change with unsed arguments - if(!PanelMode) depth += stride + offset; -#endif conj_if::IsComplex && Conjugate> cj; const_blas_data_mapper lhs(_lhs,lhsStride); Index count = 0; @@ -1219,11 +1217,9 @@ EIGEN_DONT_INLINE void gemm_pack_rhs=depth && offset<=stride)); -#ifdef __clang__ - // Workaround clang ABI change with unsed arguments - if(!PanelMode) depth += stride + offset; -#endif conj_if::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; Index count = 0; @@ -1274,11 +1270,9 @@ EIGEN_DONT_INLINE void gemm_pack_rhs=depth && offset<=stride)); -#ifdef __clang__ - // Workaround clang ABI change with unsed arguments - if(!PanelMode) depth += stride + offset; -#endif conj_if::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; Index count = 0; diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index c3e1e7142..a73ce5ff0 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h @@ -92,14 +92,8 @@ EIGEN_DONT_INLINE void general_matrix_vector_product void ignore_unused_variable(const T&) {} + } +} +#define EIGEN_UNUSED_VARIABLE(var) Eigen::internal::ignore_unused_variable(var); #if !defined(EIGEN_ASM_COMMENT) #if (defined __GNUC__) && ( defined(__i386__) || defined(__x86_64__) ) From 6b6071866bf4d7832e623ae194fd752e16832765 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Tue, 25 Feb 2014 18:55:16 +0100 Subject: [PATCH 06/24] Make pivoting HouseholderQR compatible with custom scalar types --- Eigen/src/QR/ColPivHouseholderQR.h | 2 +- Eigen/src/QR/FullPivHouseholderQR.h | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 07126a9f4..4824880f5 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -350,7 +350,7 @@ template class ColPivHouseholderQR return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. - : NumTraits::epsilon() * m_qr.diagonalSize(); + : NumTraits::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 44eaa1b1a..a7b0fc16f 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -346,7 +346,7 @@ template class FullPivHouseholderQR return m_usePrescribedThreshold ? m_prescribedThreshold // this formula comes from experimenting (see "LU precision tuning" thread on the list) // and turns out to be identical to Higham's formula used already in LDLt. - : NumTraits::epsilon() * m_qr.diagonalSize(); + : NumTraits::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. @@ -417,7 +417,7 @@ FullPivHouseholderQR& FullPivHouseholderQR::compute(cons m_temp.resize(cols); - m_precision = NumTraits::epsilon() * size; + m_precision = NumTraits::epsilon() * RealScalar(size); m_rows_transpositions.resize(size); m_cols_transpositions.resize(size); From ac69d8769f457366581b50ae8fa620634cd9aaa1 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 26 Feb 2014 10:12:27 +0100 Subject: [PATCH 07/24] Remove early termination in LDLT: the zero on the diagonal of the input matrix does not mean the matrix is not full rank. Typical examples are matrices coming from LS with linear equality constraints. --- Eigen/src/Cholesky/LDLT.h | 14 ++++---------- test/cholesky.cpp | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+), 10 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index f34d26465..b43e85e7f 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -291,13 +291,6 @@ template<> struct ldlt_inplace cutoff = abs(NumTraits::epsilon() * biggest_in_corner); } - // Finish early if the matrix is not full rank. - if(biggest_in_corner < cutoff) - { - for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; - break; - } - transpositions.coeffRef(k) = index_of_biggest_in_corner; if(k != index_of_biggest_in_corner) { @@ -333,6 +326,7 @@ template<> struct ldlt_inplace if(rs>0) A21.noalias() -= A20 * temp.head(k); } + if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) A21 /= mat.coeffRef(k,k); @@ -518,12 +512,12 @@ struct solve_retval, Rhs> typedef typename LDLTType::RealScalar RealScalar; const Diagonal vectorD = dec().vectorD(); RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits::epsilon(), - RealScalar(1) / NumTraits::highest()); // motivated by LAPACK's xGELSS + RealScalar(1) / NumTraits::highest()); // motivated by LAPACK's xGELSS for (Index i = 0; i < vectorD.size(); ++i) { if(abs(vectorD(i)) > tolerance) - dst.row(i) /= vectorD(i); + dst.row(i) /= vectorD(i); else - dst.row(i).setZero(); + dst.row(i).setZero(); } // dst = L^-T (D^-1 L^-1 P b) diff --git a/test/cholesky.cpp b/test/cholesky.cpp index b980dc572..d4d90e467 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -179,6 +179,38 @@ template void cholesky(const MatrixType& m) // restore if(sign == -1) symm = -symm; + + // check matrices coming from linear constraints with Lagrange multipliers + if(rows>=3) + { + SquareMatrixType A = symm; + int c = internal::random(0,rows-2); + A.bottomRightCorner(c,c).setZero(); + // Make sure a solution exists: + vecX.setRandom(); + vecB = A * vecX; + vecX.setZero(); + ldltlo.compute(A); + VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); + VERIFY_IS_APPROX(A * vecX, vecB); + } + + // check non-full rank matrices + if(rows>=3) + { + int r = internal::random(1,rows-1); + Matrix a = Matrix::Random(rows,r); + SquareMatrixType A = a * a.adjoint(); + // Make sure a solution exists: + vecX.setRandom(); + vecB = A * vecX; + vecX.setZero(); + ldltlo.compute(A); + VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); + VERIFY_IS_APPROX(A * vecX, vecB); + } } // update/downdate From 41e89c73c755fe70f8babea4a8a3b491057e5c39 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Thu, 27 Feb 2014 12:57:24 +0100 Subject: [PATCH 08/24] Regression test for bug #752 --- test/block.cpp | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/test/block.cpp b/test/block.cpp index 189fa5aba..9ed5d7bc5 100644 --- a/test/block.cpp +++ b/test/block.cpp @@ -10,6 +10,26 @@ #define EIGEN_NO_STATIC_ASSERT // otherwise we fail at compile time on unused paths #include "main.h" +template +typename Eigen::internal::enable_if::IsComplex,typename MatrixType::Scalar>::type +block_real_only(const MatrixType &m1, Index r1, Index r2, Index c1, Index c2, const Scalar& s1) { + // check cwise-Functions: + VERIFY_IS_APPROX(m1.row(r1).cwiseMax(s1), m1.cwiseMax(s1).row(r1)); + VERIFY_IS_APPROX(m1.col(c1).cwiseMin(s1), m1.cwiseMin(s1).col(c1)); + + VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMin(s1), m1.cwiseMin(s1).block(r1,c1,r2-r1+1,c2-c1+1)); + VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMax(s1), m1.cwiseMax(s1).block(r1,c1,r2-r1+1,c2-c1+1)); + + return Scalar(0); +} + +template +typename Eigen::internal::enable_if::IsComplex,typename MatrixType::Scalar>::type +block_real_only(const MatrixType &, Index, Index, Index, Index, const Scalar&) { + return Scalar(0); +} + + template void block(const MatrixType& m) { typedef typename MatrixType::Index Index; @@ -37,6 +57,8 @@ template void block(const MatrixType& m) Index c1 = internal::random(0,cols-1); Index c2 = internal::random(c1,cols-1); + block_real_only(m1, r1, r2, c1, c1, s1); + //check row() and col() VERIFY_IS_EQUAL(m1.col(c1).transpose(), m1.transpose().row(c1)); //check operator(), both constant and non-constant, on row() and col() @@ -51,7 +73,8 @@ template void block(const MatrixType& m) VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + s1 * m1_copy.col(c2)); m1.col(c1).col(0) += s1 * m1_copy.col(c2); VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + Scalar(2) * s1 * m1_copy.col(c2)); - + + //check block() Matrix b1(1,1); b1(0,0) = m1(r1,c1); From 76d2ca27e51179b67a246a51d9f0396991606ac8 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 28 Feb 2014 13:11:39 +0100 Subject: [PATCH 09/24] Fix PaStiX support for Pastix 5.2 --- Eigen/src/PaStiXSupport/PaStiXSupport.h | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Eigen/src/PaStiXSupport/PaStiXSupport.h b/Eigen/src/PaStiXSupport/PaStiXSupport.h index a955287d1..8a546dc2f 100644 --- a/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -12,6 +12,14 @@ namespace Eigen { +#if defined(DCOMPLEX) + #define PASTIX_COMPLEX COMPLEX + #define PASTIX_DCOMPLEX DCOMPLEX +#else + #define PASTIX_COMPLEX std::complex + #define PASTIX_DCOMPLEX std::complex +#endif + /** \ingroup PaStiXSupport_Module * \brief Interface to the PaStix solver * @@ -74,14 +82,14 @@ namespace internal { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } if (nbrhs == 0) {x = NULL; nbrhs=1;} - c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast(vals), perm, invp, reinterpret_cast(x), nbrhs, iparm, dparm); + c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast(vals), perm, invp, reinterpret_cast(x), nbrhs, iparm, dparm); } void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex *vals, int *perm, int * invp, std::complex *x, int nbrhs, int *iparm, double *dparm) { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } if (nbrhs == 0) {x = NULL; nbrhs=1;} - z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast(vals), perm, invp, reinterpret_cast(x), nbrhs, iparm, dparm); + z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast(vals), perm, invp, reinterpret_cast(x), nbrhs, iparm, dparm); } // Convert the matrix to Fortran-style Numbering From 47679c50ae2f51edaa18e4299ae361cfa774be96 Mon Sep 17 00:00:00 2001 From: Olivier Saut Date: Mon, 3 Mar 2014 14:44:39 +0100 Subject: [PATCH 10/24] Typo in the example for Eigen::SelfAdjointEigenSolver::eigenvectors, the first eigenvector should be col(0) not col(1) --- doc/snippets/EigenSolver_eigenvectors.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/snippets/EigenSolver_eigenvectors.cpp b/doc/snippets/EigenSolver_eigenvectors.cpp index 0fad4dadb..8355f76c9 100644 --- a/doc/snippets/EigenSolver_eigenvectors.cpp +++ b/doc/snippets/EigenSolver_eigenvectors.cpp @@ -1,4 +1,4 @@ MatrixXd ones = MatrixXd::Ones(3,3); EigenSolver es(ones); -cout << "The first eigenvector of the 3x3 matrix of ones is:" - << endl << es.eigenvectors().col(1) << endl; +cout << "The first eigenvector of the 3x3 matrix of ones is:" + << endl << es.eigenvectors().col(0) << endl; From 04e1e38eedad891862083de85b965b8660dea6f7 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Tue, 4 Mar 2014 15:10:29 +0100 Subject: [PATCH 11/24] bug #289: Removed useless static keywords --- Eigen/src/Core/util/Memory.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index bc472d720..7bd3056d5 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -464,7 +464,7 @@ template inline void conditional_aligned_delete_auto(T * * There is also the variant first_aligned(const MatrixBase&) defined in DenseCoeffsBase.h. */ template -static inline Index first_aligned(const Scalar* array, Index size) +inline Index first_aligned(const Scalar* array, Index size) { enum { PacketSize = packet_traits::size, PacketAlignedMask = PacketSize-1 @@ -492,7 +492,7 @@ static inline Index first_aligned(const Scalar* array, Index size) /** \internal Returns the smallest integer multiple of \a base and greater or equal to \a size */ template -inline static Index first_multiple(Index size, Index base) +inline Index first_multiple(Index size, Index base) { return ((size+base-1)/base)*base; } From 7313f32efa5fe881469e14f81268a98c6d488ee7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 4 Mar 2014 17:24:00 +0100 Subject: [PATCH 12/24] Help MSVC to inline some trivial functions --- Eigen/src/Core/ArrayWrapper.h | 2 +- Eigen/src/Core/MatrixBase.h | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Core/ArrayWrapper.h b/Eigen/src/Core/ArrayWrapper.h index 21830745b..4bb648024 100644 --- a/Eigen/src/Core/ArrayWrapper.h +++ b/Eigen/src/Core/ArrayWrapper.h @@ -49,7 +49,7 @@ class ArrayWrapper : public ArrayBase > typedef typename internal::nested::type NestedExpressionType; EIGEN_DEVICE_FUNC - inline ArrayWrapper(ExpressionType& matrix) : m_expression(matrix) {} + EIGEN_STRONG_INLINE ArrayWrapper(ExpressionType& matrix) : m_expression(matrix) {} EIGEN_DEVICE_FUNC inline Index rows() const { return m_expression.rows(); } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 49b69f873..172929562 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -356,8 +356,6 @@ template class MatrixBase Scalar trace() const; -/////////// Array module /////////// - template EIGEN_DEVICE_FUNC RealScalar lpNorm() const; EIGEN_DEVICE_FUNC MatrixBase& matrix() { return *this; } @@ -365,8 +363,10 @@ template class MatrixBase /** \returns an \link Eigen::ArrayBase Array \endlink expression of this matrix * \sa ArrayBase::matrix() */ - EIGEN_DEVICE_FUNC ArrayWrapper array() { return derived(); } - EIGEN_DEVICE_FUNC const ArrayWrapper array() const { return derived(); } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ArrayWrapper array() { return derived(); } + /** \returns a const \link Eigen::ArrayBase Array \endlink expression of this matrix + * \sa ArrayBase::matrix() */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArrayWrapper array() const { return derived(); } /////////// LU module /////////// From d5cc083782bc59049755696f2d22edfb0172d124 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Wed, 5 Mar 2014 14:50:00 +0100 Subject: [PATCH 13/24] Fixed bug #754. Only inserted (!defined(_WIN32_WCE)) analog to alloc and free implementation (not tested, but should be correct). --- Eigen/src/Core/util/Memory.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 7bd3056d5..286e963d2 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -274,12 +274,12 @@ inline void* aligned_realloc(void *ptr, size_t new_size, size_t old_size) // The defined(_mm_free) is just here to verify that this MSVC version // implements _mm_malloc/_mm_free based on the corresponding _aligned_ // functions. This may not always be the case and we just try to be safe. - #if defined(_MSC_VER) && defined(_mm_free) + #if defined(_MSC_VER) && (!defined(_WIN32_WCE)) && defined(_mm_free) result = _aligned_realloc(ptr,new_size,16); #else result = generic_aligned_realloc(ptr,new_size,old_size); #endif -#elif defined(_MSC_VER) +#elif defined(_MSC_VER) && (!defined(_WIN32_WCE)) result = _aligned_realloc(ptr,new_size,16); #else result = handmade_aligned_realloc(ptr,new_size,old_size); From ce41b72eb8b0657957893275e4936b0961bae80b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 7 Mar 2014 23:09:39 +0100 Subject: [PATCH 14/24] Extend sizeof unit test --- test/sizeof.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/sizeof.cpp b/test/sizeof.cpp index d9ad35620..7044d2062 100644 --- a/test/sizeof.cpp +++ b/test/sizeof.cpp @@ -21,6 +21,8 @@ template void verifySizeOf(const MatrixType&) void test_sizeof() { CALL_SUBTEST(verifySizeOf(Matrix()) ); + CALL_SUBTEST(verifySizeOf(Vector2d()) ); + CALL_SUBTEST(verifySizeOf(Vector4f()) ); CALL_SUBTEST(verifySizeOf(Matrix4d()) ); CALL_SUBTEST(verifySizeOf(Matrix()) ); CALL_SUBTEST(verifySizeOf(Matrix()) ); From 33ca9b4ee690ed12bd072273ea0d24c5eaae77fd Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 7 Mar 2014 23:11:38 +0100 Subject: [PATCH 15/24] Add support for OSX in BTL and fix a few warnings --- bench/btl/CMakeLists.txt | 8 +++- bench/btl/generic_bench/bench.hh | 4 +- bench/btl/generic_bench/btl.hh | 2 +- bench/btl/generic_bench/init/init_matrix.hh | 10 ++-- bench/btl/generic_bench/init/init_vector.hh | 2 +- .../timers/portable_perf_analyzer.hh | 2 +- .../generic_bench/timers/portable_timer.hh | 46 ++++++++++++++++++- bench/btl/libs/eigen3/eigen3_interface.hh | 8 ++-- 8 files changed, 65 insertions(+), 17 deletions(-) diff --git a/bench/btl/CMakeLists.txt b/bench/btl/CMakeLists.txt index 119b470d9..0ac5c15fc 100644 --- a/bench/btl/CMakeLists.txt +++ b/bench/btl/CMakeLists.txt @@ -48,6 +48,12 @@ include_directories( # set(DEFAULT_LIBRARIES ${MKL_LIBRARIES}) # endif (MKL_FOUND) +find_library(EIGEN_BTL_RT_LIBRARY rt) +# if we cannot find it easily, then we don't need it! +if(NOT EIGEN_BTL_RT_LIBRARY) + set(EIGEN_BTL_RT_LIBRARY "") +endif() + MACRO(BTL_ADD_BENCH targetname) foreach(_current_var ${ARGN}) @@ -70,7 +76,7 @@ MACRO(BTL_ADD_BENCH targetname) IF(BUILD_${targetname}) ADD_EXECUTABLE(${targetname} ${_sources}) ADD_TEST(${targetname} "${targetname}") - target_link_libraries(${targetname} ${DEFAULT_LIBRARIES} rt) + target_link_libraries(${targetname} ${DEFAULT_LIBRARIES} ${EIGEN_BTL_RT_LIBRARY}) ENDIF(BUILD_${targetname}) ENDMACRO(BTL_ADD_BENCH) diff --git a/bench/btl/generic_bench/bench.hh b/bench/btl/generic_bench/bench.hh index 005c36395..7b7b951b5 100644 --- a/bench/btl/generic_bench/bench.hh +++ b/bench/btl/generic_bench/bench.hh @@ -102,8 +102,8 @@ BTL_DONT_INLINE void bench( int size_min, int size_max, int nb_point ) // merge the two data std::vector newSizes; std::vector newFlops; - int i=0; - int j=0; + unsigned int i=0; + unsigned int j=0; while (i BTL_DONT_INLINE void init_matrix(Vector & A, int size){ A.resize(size); - for (int row=0; row(A[row],size,row); } } @@ -50,11 +50,11 @@ BTL_DONT_INLINE void init_matrix(Vector & A, int size){ template BTL_DONT_INLINE void init_matrix_symm(Matrix& A, int size){ A.resize(size); - for (int row=0; row @@ -87,6 +87,48 @@ }; // Portable_Timer +#elif defined(__APPLE__) +#include +#include + + +class Portable_Timer +{ + public: + + Portable_Timer() + { + } + + void start() + { + m_start_time = double(mach_absolute_time())*1e-9;; + + } + + void stop() + { + m_stop_time = double(mach_absolute_time())*1e-9;; + + } + + double elapsed() + { + return user_time(); + } + + double user_time() + { + return m_stop_time - m_start_time; + } + + +private: + + double m_stop_time, m_start_time; + +}; // Portable_Timer (Apple) + #else #include @@ -138,7 +180,7 @@ private: int m_clkid; double m_stop_time, m_start_time; -}; // Portable_Timer +}; // Portable_Timer (Linux) #endif diff --git a/bench/btl/libs/eigen3/eigen3_interface.hh b/bench/btl/libs/eigen3/eigen3_interface.hh index 31bcc1f93..0c8aa74da 100644 --- a/bench/btl/libs/eigen3/eigen3_interface.hh +++ b/bench/btl/libs/eigen3/eigen3_interface.hh @@ -52,8 +52,8 @@ public : static BTL_DONT_INLINE void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ A.resize(A_stl[0].size(), A_stl.size()); - for (int j=0; j Date: Fri, 7 Mar 2014 23:13:14 +0100 Subject: [PATCH 16/24] Fix typo and formating --- Eigen/src/Geometry/Scaling.h | 8 ++++---- doc/TopicMultithreading.dox | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Geometry/Scaling.h b/Eigen/src/Geometry/Scaling.h index 1c25f36fe..023fba2ee 100644 --- a/Eigen/src/Geometry/Scaling.h +++ b/Eigen/src/Geometry/Scaling.h @@ -62,10 +62,10 @@ public: template inline Transform operator* (const Transform& t) const { - Transform res = t; - res.prescale(factor()); - return res; -} + Transform res = t; + res.prescale(factor()); + return res; + } /** Concatenates a uniform scaling and a linear transformation matrix */ // TODO returns an expression diff --git a/doc/TopicMultithreading.dox b/doc/TopicMultithreading.dox index fb944af29..ba5e26290 100644 --- a/doc/TopicMultithreading.dox +++ b/doc/TopicMultithreading.dox @@ -39,7 +39,7 @@ int main(int argc, char** argv) } \endcode -\warning note that all functions generating random matrices are \b not re-entrant nor thread-safe. Those include DenseBase::Random(), and DenseBase::setRandom() despite a call to Eigen::initParallel(). This is because these functions are based on std::rand which is not re-entrant. For thread-safe random generator, we recommend the use of boost::random of c++11 random feature. +\warning note that all functions generating random matrices are \b not re-entrant nor thread-safe. Those include DenseBase::Random(), and DenseBase::setRandom() despite a call to Eigen::initParallel(). This is because these functions are based on std::rand which is not re-entrant. For thread-safe random generator, we recommend the use of boost::random or c++11 random feature. In the case your application is parallelized with OpenMP, you might want to disable Eigen's own parallization as detailed in the previous section. From ba2f79e6800ca2762821fe5ab537bc5f10a1268a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 7 Mar 2014 23:18:20 +0100 Subject: [PATCH 17/24] Fix selfadjoint_matrix_vector_product for complex with packet size > 2 (e.g., AVX) --- Eigen/src/Core/products/SelfadjointMatrixVector.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index f698f67f9..fdc81205a 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -113,9 +113,9 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product huge speed up) // gcc 4.2 does this optimization automatically. From 033ee7f6d91286920d56de1003f582f4570ff5b3 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Sat, 8 Mar 2014 00:44:56 -0500 Subject: [PATCH 18/24] Fix typo: 'explicitely' -> 'explicitly' --- doc/TopicLinearAlgebraDecompositions.dox | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/TopicLinearAlgebraDecompositions.dox b/doc/TopicLinearAlgebraDecompositions.dox index 8649cc27b..77f2c92ab 100644 --- a/doc/TopicLinearAlgebraDecompositions.dox +++ b/doc/TopicLinearAlgebraDecompositions.dox @@ -249,7 +249,7 @@ For an introduction on linear solvers and decompositions, check this \link Tutor
Implicit Multi Threading (MT)
Means the algorithm can take advantage of multicore processors via OpenMP. "Implicit" means the algortihm itself is not parallelized, but that it relies on parallelized matrix-matrix product rountines.
Explicit Multi Threading (MT)
-
Means the algorithm is explicitely parallelized to take advantage of multicore processors via OpenMP.
+
Means the algorithm is explicitly parallelized to take advantage of multicore processors via OpenMP.
Meta-unroller
Means the algorithm is automatically and explicitly unrolled for very small fixed size matrices.
From ce99b502ce1d80129a028ddf02fff51f6c51249b Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Tue, 17 Dec 2013 10:49:43 -0800 Subject: [PATCH 19/24] Use vectorization when packing row-major rhs matrices. (bug #717) --- .../Core/products/GeneralBlockPanelKernel.h | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 08cc14bd7..686ff84f1 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -1261,6 +1261,7 @@ EIGEN_DONT_INLINE void gemm_pack_rhs struct gemm_pack_rhs { + typedef typename packet_traits::type Packet; enum { PacketSize = packet_traits::size }; EIGEN_DONT_INLINE void operator()(Scalar* blockB, const Scalar* rhs, Index rhsStride, Index depth, Index cols, Index stride=0, Index offset=0); }; @@ -1282,12 +1283,18 @@ EIGEN_DONT_INLINE void gemm_pack_rhs(&rhs[k*rhsStride + j2]); + pstoreu(blockB+count, cj.pconj(A)); + count += PacketSize; + } else { + const Scalar* b0 = &rhs[k*rhsStride + j2]; + blockB[count+0] = cj(b0[0]); + blockB[count+1] = cj(b0[1]); + if(nr==4) blockB[count+2] = cj(b0[2]); + if(nr==4) blockB[count+3] = cj(b0[3]); + count += nr; + } } // skip what we have after if(PanelMode) count += nr * (stride-offset-depth); From a6b130c63c4e1cb6be6bf71414f80feaea922b23 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 5 Mar 2014 10:07:44 +0100 Subject: [PATCH 20/24] swap 3.2 <-> default CTestConfig.cmake file --- CTestConfig.cmake | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/CTestConfig.cmake b/CTestConfig.cmake index 0557c491a..4c0027824 100644 --- a/CTestConfig.cmake +++ b/CTestConfig.cmake @@ -4,10 +4,14 @@ ## # The following are required to uses Dart and the Cdash dashboard ## ENABLE_TESTING() ## INCLUDE(CTest) -set(CTEST_PROJECT_NAME "Eigen3.2") +set(CTEST_PROJECT_NAME "Eigen") set(CTEST_NIGHTLY_START_TIME "00:00:00 UTC") set(CTEST_DROP_METHOD "http") set(CTEST_DROP_SITE "manao.inria.fr") -set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen3.2") +set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen") set(CTEST_DROP_SITE_CDASH TRUE) +set(CTEST_PROJECT_SUBPROJECTS +Official +Unsupported +) From 804ef2350d51f6fe4d22f18036dfc7e10a4bb711 Mon Sep 17 00:00:00 2001 From: Abraham Bachrach Date: Sun, 9 Mar 2014 16:56:44 -0700 Subject: [PATCH 21/24] Move the Base typedef's from private to public scope Move the Quaternion::Base typedef from private to public scope so that one may create child classes of Quaternion. NOTE: This matches the semantics of MatrixBase. --- Eigen/src/Geometry/Quaternion.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 803001272..6ee0e4f75 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -34,8 +34,9 @@ struct quaternionbase_assign_impl; template class QuaternionBase : public RotationBase { + public: typedef RotationBase Base; -public: + using Base::operator*; using Base::derived; @@ -223,10 +224,10 @@ struct traits > template class Quaternion : public QuaternionBase > { +public: typedef QuaternionBase > Base; enum { IsAligned = internal::traits::IsAligned }; -public: typedef _Scalar Scalar; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Quaternion) @@ -334,9 +335,9 @@ template class Map, _Options > : public QuaternionBase, _Options> > { + public: typedef QuaternionBase, _Options> > Base; - public: typedef _Scalar Scalar; typedef typename internal::traits::Coefficients Coefficients; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Map) @@ -371,9 +372,9 @@ template class Map, _Options > : public QuaternionBase, _Options> > { + public: typedef QuaternionBase, _Options> > Base; - public: typedef _Scalar Scalar; typedef typename internal::traits::Coefficients Coefficients; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Map) From bbc0ada12a3ae1531860205b7f1e21f991bb6273 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Tue, 11 Mar 2014 12:18:32 +0100 Subject: [PATCH 22/24] Avoid stupid "enumeral mismatch in conditional expression" warnings in GCC --- Eigen/src/Geometry/Transform.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 0a5012cfe..f40644011 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -530,9 +530,9 @@ public: inline Transform& operator=(const UniformScaling& t); inline Transform& operator*=(const UniformScaling& s) { return scale(s.factor()); } - inline Transform operator*(const UniformScaling& s) const + inline TransformTimeDiagonalReturnType operator*(const UniformScaling& s) const { - Transform res = *this; + TransformTimeDiagonalReturnType res = *this; res.scale(s.factor()); return res; } From 88aa18df641f8235c88661e41ea128f33d88a8f6 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Wed, 12 Mar 2014 13:43:19 +0100 Subject: [PATCH 23/24] bug #759: Removed hard-coded double-math from Quaternion::angularDistance. Some documentation improvements --- Eigen/src/Geometry/Quaternion.h | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 6ee0e4f75..a712692e5 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -204,6 +204,8 @@ class QuaternionBase : public RotationBase * \li \c Quaternionf for \c float * \li \c Quaterniond for \c double * + * \warning Operations interpreting the quaternion as rotation have undefined behavior if the quaternion is not normalized. + * * \sa class AngleAxis, class Transform */ @@ -345,7 +347,7 @@ class Map, _Options > /** Constructs a Mapped Quaternion object from the pointer \a coeffs * - * The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order: + * The pointer \a coeffs must reference the four coefficients of Quaternion in the following order: * \code *coeffs == {x, y, z, w} \endcode * * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */ @@ -465,7 +467,7 @@ QuaternionBase::_transformVector(Vector3 v) const // Note that this algorithm comes from the optimization by hand // of the conversion to a Matrix followed by a Matrix/Vector product. // It appears to be much faster than the common algorithm found - // in the litterature (30 versus 39 flops). It also requires two + // in the literature (30 versus 39 flops). It also requires two // Vector3 as temporaries. Vector3 uv = this->vec().cross(v); uv += uv; @@ -668,10 +670,10 @@ QuaternionBase::angularDistance(const QuaternionBase& oth { using std::acos; using std::abs; - double d = abs(this->dot(other)); - if (d>=1.0) + Scalar d = abs(this->dot(other)); + if (d>=Scalar(1)) return Scalar(0); - return static_cast(2 * acos(d)); + return Scalar(2) * acos(d); } From 2379ccffcb8f1af10f8ed91353d5cb8b3a9a7847 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Wed, 12 Mar 2014 13:48:09 +0100 Subject: [PATCH 24/24] bug #755: CommaInitializer produced wrong assertions in absence of ReturnValueOptimization. --- Eigen/src/Core/CommaInitializer.h | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/CommaInitializer.h b/Eigen/src/Core/CommaInitializer.h index 2bbf74b05..70cbfeff5 100644 --- a/Eigen/src/Core/CommaInitializer.h +++ b/Eigen/src/Core/CommaInitializer.h @@ -45,6 +45,18 @@ struct CommaInitializer m_xpr.block(0, 0, other.rows(), other.cols()) = other; } + /* Copy/Move constructor which transfers ownership. This is crucial in + * absence of return value optimization to avoid assertions during destruction. */ + // FIXME in C++11 mode this could be replaced by a proper RValue constructor + EIGEN_DEVICE_FUNC + inline CommaInitializer(const CommaInitializer& o) + : m_xpr(o.m_xpr), m_row(o.m_row), m_col(o.m_col), m_currentBlockRows(o.m_currentBlockRows) { + // Mark original object as finished. In absence of R-value references we need to const_cast: + const_cast(o).m_row = m_xpr.rows(); + const_cast(o).m_col = m_xpr.cols(); + const_cast(o).m_currentBlockRows = 0; + } + /* inserts a scalar value in the target matrix */ EIGEN_DEVICE_FUNC CommaInitializer& operator,(const Scalar& s) @@ -110,7 +122,7 @@ struct CommaInitializer EIGEN_DEVICE_FUNC inline XprType& finished() { return m_xpr; } - XprType& m_xpr; // target expression + XprType& m_xpr; // target expression Index m_row; // current row id Index m_col; // current col id Index m_currentBlockRows; // current block height