From 544888e342bb7d75412bf0e17b66a3c6196d24eb Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 7 Jul 2009 09:05:20 +0200 Subject: [PATCH] add a generic mechanism to copy a special matrix to a dense matrix so that we don't need to add other specialization of MatrixBase::operator=, Matrix::=, and Matrix::Matrix(...) --- Eigen/src/Cholesky/LLT.h | 41 ----------------------------- Eigen/src/Core/DiagonalMatrix.h | 11 ++++---- Eigen/src/Core/Matrix.h | 30 ++++++--------------- Eigen/src/Core/MatrixBase.h | 24 ++++++++--------- Eigen/src/Core/TriangularMatrix.h | 33 ++++++++++++----------- Eigen/src/Sparse/SparseMatrixBase.h | 28 ++++++++------------ 6 files changed, 54 insertions(+), 113 deletions(-) diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index bad360579..fccd53459 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -225,47 +225,6 @@ template struct LLT_Traits template void LLT::compute(const MatrixType& a) { -// assert(a.rows()==a.cols()); -// const int size = a.rows(); -// m_matrix.resize(size, size); -// // The biggest overall is the point of reference to which further diagonals -// // are compared; if any diagonal is negligible compared -// // to the largest overall, the algorithm bails. This cutoff is suggested -// // 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. -// const RealScalar cutoff = machine_epsilon() * size * a.diagonal().cwise().abs().maxCoeff(); -// RealScalar x; -// x = ei_real(a.coeff(0,0)); -// m_matrix.coeffRef(0,0) = ei_sqrt(x); -// if(size==1) -// { -// m_isInitialized = true; -// return; -// } -// m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0)); -// for (int j = 1; j < size; ++j) -// { -// x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); -// if (ei_abs(x) < cutoff) continue; -// -// m_matrix.coeffRef(j,j) = x = ei_sqrt(x); -// -// int endSize = size-j-1; -// if (endSize>0) { -// // Note that when all matrix columns have good alignment, then the following -// // product is guaranteed to be optimal with respect to alignment. -// m_matrix.col(j).end(endSize) = -// (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy(); -// -// // FIXME could use a.col instead of a.row -// m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint() -// - m_matrix.col(j).end(endSize) ) / x; -// } -// } -// -// m_isInitialized = true; - assert(a.rows()==a.cols()); const int size = a.rows(); m_matrix.resize(size, size); diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h index a092b7794..5fc80c92b 100644 --- a/Eigen/src/Core/DiagonalMatrix.h +++ b/Eigen/src/Core/DiagonalMatrix.h @@ -50,6 +50,8 @@ class DiagonalBase : public MultiplierBase #endif // not EIGEN_PARSED_BY_DOXYGEN DenseMatrixType toDenseMatrix() const { return derived(); } + template + void evalToDense(MatrixBase &other) const; inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); } inline DiagonalVectorType& diagonal() { return derived().diagonal(); } @@ -63,12 +65,11 @@ class DiagonalBase : public MultiplierBase }; template -template -Derived& MatrixBase::operator=(const DiagonalBase &other) +template +void DiagonalBase::evalToDense(MatrixBase &other) const { - setZero(); - diagonal() = other.diagonal(); - return derived(); + other.setZero(); + other.diagonal() = diagonal(); } /** \class DiagonalMatrix diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index b7a8ebf21..18c8f969a 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -439,34 +439,20 @@ class Matrix { other.evalTo(*this); } /** Destructor */ inline ~Matrix() {} - - template - EIGEN_STRONG_INLINE Matrix& operator=(const DiagonalBase &other) + /** \sa MatrixBase::operator=(const AnyMatrixBase&) */ + template + EIGEN_STRONG_INLINE Matrix& operator=(const AnyMatrixBase &other) { - resize(other.diagonal().size(), other.diagonal().size()); - Base::operator=(other); - return *this; - } - - template - EIGEN_STRONG_INLINE Matrix(const DiagonalBase &other) - : m_storage(other.diagonal().size() * other.diagonal().size(), other.diagonal().size(), other.diagonal().size()) - { - *this = other; - } - - template - EIGEN_STRONG_INLINE Matrix& operator=(const TriangularBase &other) - { - resize(other.rows(), other.cols()); + resize(other.derived().rows(), other.derived().cols()); Base::operator=(other.derived()); return *this; } - template - EIGEN_STRONG_INLINE Matrix(const TriangularBase &other) - : m_storage(other.rows() * other.cols(), other.rows(), other.cols()) + /** \sa MatrixBase::operator=(const AnyMatrixBase&) */ + template + EIGEN_STRONG_INLINE Matrix(const AnyMatrixBase &other) + : m_storage(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols()) { *this = other; } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 82417a144..7aedbd3c7 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -284,6 +284,17 @@ template class MatrixBase */ Derived& operator=(const MatrixBase& other); + /** Copies the generic expression \a other into *this. \returns a reference to *this. + * The expression must provide a (templated) evalToDense(Derived& dst) const function + * which does the actual job. In practice, this allows any user to write its own + * special matrix without having to modify MatrixBase */ + template + Derived& operator=(const AnyMatrixBase &other) + { other.derived().evalToDense(derived()); return derived(); } + + template + Derived& operator=(const ReturnByValue& func); + #ifndef EIGEN_PARSED_BY_DOXYGEN /** Copies \a other into *this without evaluating other. \returns a reference to *this. */ template @@ -297,10 +308,6 @@ template class MatrixBase template Derived& lazyAssign(const Flagged& other) { return lazyAssign(other._expression()); } - - /** Overloaded for fast triangular part to dense matrix evaluation */ - template - Derived& lazyAssign(const TriangularBase &other); #endif // not EIGEN_PARSED_BY_DOXYGEN CommaInitializer operator<< (const Scalar& s); @@ -403,12 +410,6 @@ template class MatrixBase const DiagonalProduct operator*(const DiagonalBase &diagonal) const; - template - Derived& operator=(const DiagonalBase &other); - - template - Derived& operator=(const TriangularBase &other); - template typename ei_plain_matrix_type_column_major::type solveTriangular(const MatrixBase& other) const; @@ -738,9 +739,6 @@ template class MatrixBase // dense = dense * sparse template Derived& lazyAssign(const SparseProduct& product); - - template - Derived& operator=(const ReturnByValue& func); #ifdef EIGEN_MATRIXBASE_PLUGIN #include EIGEN_MATRIXBASE_PLUGIN diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 8131ef323..81621dbcc 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -90,6 +90,11 @@ template class TriangularBase : public MultiplierBase inline Derived& derived() { return *static_cast(this); } #endif // not EIGEN_PARSED_BY_DOXYGEN + template + void evalToDense(MatrixBase &other) const; + template + void evalToDenseLazy(MatrixBase &other) const; + protected: void check_coordinates(int row, int col) @@ -516,36 +521,34 @@ void TriangularView::lazyAssign(const TriangularBase -template -Derived& MatrixBase::operator=(const TriangularBase &other) +template +void TriangularBase::evalToDense(MatrixBase &other) const { - if(ei_traits::Flags & EvalBeforeAssigningBit) + if(ei_traits::Flags & EvalBeforeAssigningBit) { - typename TriangularDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); - other_evaluated.lazyAssign(other); - this->swap(other_evaluated); + typename Derived::PlainMatrixType other_evaluated(rows(), cols()); + evalToDenseLazy(other_evaluated); + other.derived().swap(other_evaluated); } else - lazyAssign(other.derived()); - return derived(); + evalToDenseLazy(other.derived()); } /** Assigns a triangular or selfadjoint matrix to a dense matrix. * If the matrix is triangular, the opposite part is set to zero. */ template -template -Derived& MatrixBase::lazyAssign(const TriangularBase &other) +template +void TriangularBase::evalToDenseLazy(MatrixBase &other) const { - const bool unroll = Derived::SizeAtCompileTime * TriangularDerived::CoeffReadCost / 2 + const bool unroll = DenseDerived::SizeAtCompileTime * Derived::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT; ei_assert(this->rows() == other.rows() && this->cols() == other.cols()); ei_part_assignment_impl - ::ExpressionType, TriangularDerived::Mode, - unroll ? int(Derived::SizeAtCompileTime) : Dynamic, + ::ExpressionType, Derived::Mode, + unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic, true // clear the opposite triangular part - >::run(this->const_cast_derived(), other.derived()._expression()); - return derived(); + >::run(other.derived(), derived()._expression()); } /** \deprecated use MatrixBase::triangularView() */ diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 8c53757b0..82527f7ef 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -451,25 +451,19 @@ template class SparseMatrixBase // Derived& setRandom(); // Derived& setIdentity(); + template + void evalToDense(MatrixBase& dst) + { + dst.resize(rows(),cols()); + dst.setZero(); + for (int j=0; j toDense() const { - Matrix res(rows(),cols()); - res.setZero(); - for (int j=0; j