diff options
| -rw-r--r-- | Eigen/Core | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/Assign.h | 7 | ||||
| -rwxr-xr-x | Eigen/src/Core/Triangular.h | 67 | ||||
| -rw-r--r-- | Eigen/src/Core/TriangularAssign.h | 104 |
4 files changed, 120 insertions, 59 deletions
diff --git a/Eigen/Core b/Eigen/Core index e5ec8c08c..b22c0a0ba 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -52,6 +52,7 @@ namespace Eigen { #include "src/Core/Swap.h" #include "src/Core/CommaInitializer.h" #include "src/Core/Triangular.h" +#include "src/Core/TriangularAssign.h" } // namespace Eigen diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 1590fef81..2f15e2775 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -103,7 +103,8 @@ bool Vectorize = (Derived::Flags & OtherDerived::Flags & VectorizableBit) && ( (Derived::Flags & OtherDerived::Flags & Like1DArrayBit) ||((Derived::Flags&RowMajorBit) ? Derived::ColsAtCompileTime!=Dynamic && (Derived::ColsAtCompileTime%ei_packet_traits<typename Derived::Scalar>::size==0) - : Derived::RowsAtCompileTime!=Dynamic && (Derived::RowsAtCompileTime%ei_packet_traits<typename Derived::Scalar>::size==0)) )> + : Derived::RowsAtCompileTime!=Dynamic && (Derived::RowsAtCompileTime%ei_packet_traits<typename Derived::Scalar>::size==0)) ), +bool TriangularAssign = false> struct ei_assignment_impl; template<typename Derived> @@ -112,7 +113,9 @@ Derived& MatrixBase<Derived> ::lazyAssign(const MatrixBase<OtherDerived>& other) { // std::cout << "lazyAssign = " << Derived::Flags << " " << OtherDerived::Flags << "\n"; - ei_assignment_impl<Derived,OtherDerived>::execute(derived(),other.derived()); + ei_assignment_impl<Derived, OtherDerived, + Derived::Flags & (NullLowerBit | NullUpperBit)> + ::execute(derived(),other.derived()); return derived(); } diff --git a/Eigen/src/Core/Triangular.h b/Eigen/src/Core/Triangular.h index 9bc8def38..d11ac37ff 100755 --- a/Eigen/src/Core/Triangular.h +++ b/Eigen/src/Core/Triangular.h @@ -96,49 +96,6 @@ template<int Mode, typename MatrixType> class Triangular return Triangular<(Upper | Lower) xor Mode, Transpose<MatrixType> >((m_matrix.transpose())); } -#if 0 - - template<typename OtherDerived> - Triangular& operator=(const MatrixBase<OtherDerived>& other); - - /** Overloaded to provide optimal evaluation loops */ - template<typename OtherDerived> - Triangular& operator +=(const MatrixBase<OtherDerived>& other) - { - return *this = m_matrix + other; - } - - /** Overloaded to provide optimal evaluation loops */ - template<typename OtherDerived> - Triangular& operator *=(const MatrixBase<OtherDerived>& other) - { - return *this = this->lazyProduct(other).eval(); - } - - /** Optimized triangular matrix - matrix product */ - template<typename OtherDerived> - TriangularProduct<Mode, MatrixType, OtherDerived> lazyProduct(const MatrixBase<Scalar, OtherDerived>& other) const - { - return TriangularProduct<Mode,MatrixType,OtherDerived>(m_matrix, other.ref()); - } - - /** Optimized triangular matrix - matrix product */ - template<typename OtherDerived> - Eval<TriangularProduct<Mode, MatrixType, OtherDerived> > operator * (const MatrixBase<Scalar, OtherDerived>& other) const - { - return this->lazyProduct(other).eval(); - } - - /** Optimized matrix - triangular matrix product */ - template<typename OtherDerived> - friend Eval<Transpose<TriangularProduct<0x1 xor Mode, Transpose<MatRef>, Transpose<OtherDerived> > > > - operator * (const MatrixBase<Scalar, OtherDerived>& other, const Triangular<Mode,MatrixType>& tri) - { - return tri.transpose().lazyProduct(other.transpose()).transpose().eval(); - } - -#endif - /** \returns the product of the inverse of *this with \a other. * * This function computes the inverse-matrix matrix product inv(*this) \a other @@ -159,36 +116,32 @@ template<int Mode, typename MatrixType> class Triangular { // forward substitution if (Flags & UnitDiagBit) - res.col(c)[0] = other.col(c)[0]; + res(0,c) = other(0,c); else - res.col(c)[0] = other.col(c)[0]/_coeff(0, 0); + res(0,c) = other(0,c)/_coeff(0, 0); for (int i=1 ; i<_rows() ; ++i) { - Scalar tmp = other.col(c)[i]; - for (int j = 0 ; j < i ; ++j) - tmp -= _coeff(i,j) * res.col(c)[j]; + Scalar tmp = other(i,c) - ((this->row(i).start(i)).transpose() * res.col(c).start(i))(0,0); if (Flags & UnitDiagBit) - res.col(c)[i] = tmp; + res(i,c) = tmp; else - res.col(c)[i] = tmp/_coeff(i,i); + res(i,c) = tmp/_coeff(i,i); } } else { // backward substitution if (Flags & UnitDiagBit) - res.col(c)[_cols()-1] = other.col(c)[_cols()-1]; + res(_cols()-1,c) = other(_cols()-1,c); else - res.col(c)[_cols()-1] = other.col(c)[_cols()-1]/_coeff(_rows()-1, _cols()-1); + res(_cols()-1,c) = other(_cols()-1, c)/_coeff(_rows()-1, _cols()-1); for (int i=_rows()-2 ; i>=0 ; --i) { - Scalar tmp = other.col(c)[i]; - for (int j = i+1 ; j < _cols() ; ++j) - tmp -= _coeff(i,j) * res.col(c)[j]; + Scalar tmp = other(i,c) - ((this->row(i).end(_cols()-i-1)).transpose() * res.col(c).end(_cols()-i-1))(0,0); if (Flags & UnitDiagBit) - res.col(c)[i] = tmp; + res(i,c) = tmp; else - res.col(c)[i] = tmp/_coeff(i,i); + res(i,c) = tmp/_coeff(i,i); } } } diff --git a/Eigen/src/Core/TriangularAssign.h b/Eigen/src/Core/TriangularAssign.h new file mode 100644 index 000000000..c52085ef1 --- /dev/null +++ b/Eigen/src/Core/TriangularAssign.h @@ -0,0 +1,104 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_TRIANGULAR_ASSIGN_H +#define EIGEN_TRIANGULAR_ASSIGN_H + +template<typename Derived1, typename Derived2, int UnrollCount, int Mode> +struct ei_triangular_assign_unroller +{ + enum { + col = (UnrollCount-1) / Derived1::RowsAtCompileTime, + row = (UnrollCount-1) % Derived1::RowsAtCompileTime + }; + + static void run(Derived1 &dst, const Derived2 &src) + { + ei_triangular_assign_unroller<Derived1, Derived2, + (Mode & Lower) ? + ((row==col) ? UnrollCount-1-row : UnrollCount-1) + : ((row==0) ? UnrollCount-1-Derived1::ColsAtCompileTime+col : UnrollCount-1), + Mode>::run(dst, src); + dst.coeffRef(row, col) = src.coeff(row, col); + } +}; + +template<typename Derived1, typename Derived2, int Mode> +struct ei_triangular_assign_unroller<Derived1, Derived2, 1, Mode> +{ + static void run(Derived1 &dst, const Derived2 &src) + { + dst.coeffRef(0, 0) = src.coeff(0, 0); + } +}; + +// prevent buggy user code from causing an infinite recursion +template<typename Derived1, typename Derived2, int Mode> +struct ei_triangular_assign_unroller<Derived1, Derived2, 0, Mode> +{ + static void run(Derived1 &, const Derived2 &) {} +}; + +template<typename Derived1, typename Derived2, int Mode> +struct ei_triangular_assign_unroller<Derived1, Derived2, Dynamic, Mode> +{ + static void run(Derived1 &, const Derived2 &) {} +}; + + +template <typename Derived, typename OtherDerived, bool DummyVectorize> +struct ei_assignment_impl<Derived, OtherDerived, DummyVectorize, true> +{ + static void execute(Derived & dst, const OtherDerived & src) + { + assert(src.rows()==src.cols()); + assert(dst.rows() == src.rows() && dst.cols() == src.cols()); + + const bool unroll = Derived::SizeAtCompileTime * OtherDerived::CoeffReadCost <= EIGEN_UNROLLING_LIMIT; + + if(unroll) + { + ei_triangular_assign_unroller + <Derived, OtherDerived, unroll ? Derived::SizeAtCompileTime : Dynamic, Derived::Flags>::run + (dst.derived(), src.derived()); + } + else + { + if (Derived::Flags & Lower) + { + for(int j = 0; j < dst.cols(); j++) + for(int i = j; i < dst.rows(); i++) + dst.coeffRef(i, j) = src.coeff(i, j); + } + else + { + for(int j = 0; j < dst.cols(); j++) + for(int i = 0; i <= j; i++) + dst.coeffRef(i, j) = src.coeff(i, j); + } + } + } +}; + +#endif // EIGEN_TRIANGULAR_ASSIGN_H |