From b4c974d0597546e255edb3758a0bc58515674692 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 26 Apr 2008 19:20:26 +0000 Subject: [PATCH] Added triangular assignement, e.g.: m.upper() = a+b; only updates the upper triangular part of m. Note that: m = (a+b).upper(); updates all coefficients of m (but half of the additions will be skiped) Updated back/forward substitution to better use Eigen's capability. --- Eigen/Core | 1 + Eigen/src/Core/Assign.h | 7 +- Eigen/src/Core/Triangular.h | 67 +++---------------- Eigen/src/Core/TriangularAssign.h | 104 ++++++++++++++++++++++++++++++ 4 files changed, 120 insertions(+), 59 deletions(-) create mode 100644 Eigen/src/Core/TriangularAssign.h 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::size==0) - : Derived::RowsAtCompileTime!=Dynamic && (Derived::RowsAtCompileTime%ei_packet_traits::size==0)) )> + : Derived::RowsAtCompileTime!=Dynamic && (Derived::RowsAtCompileTime%ei_packet_traits::size==0)) ), +bool TriangularAssign = false> struct ei_assignment_impl; template @@ -112,7 +113,9 @@ Derived& MatrixBase ::lazyAssign(const MatrixBase& other) { // std::cout << "lazyAssign = " << Derived::Flags << " " << OtherDerived::Flags << "\n"; - ei_assignment_impl::execute(derived(),other.derived()); + ei_assignment_impl + ::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 class Triangular return Triangular<(Upper | Lower) xor Mode, Transpose >((m_matrix.transpose())); } -#if 0 - - template - Triangular& operator=(const MatrixBase& other); - - /** Overloaded to provide optimal evaluation loops */ - template - Triangular& operator +=(const MatrixBase& other) - { - return *this = m_matrix + other; - } - - /** Overloaded to provide optimal evaluation loops */ - template - Triangular& operator *=(const MatrixBase& other) - { - return *this = this->lazyProduct(other).eval(); - } - - /** Optimized triangular matrix - matrix product */ - template - TriangularProduct lazyProduct(const MatrixBase& other) const - { - return TriangularProduct(m_matrix, other.ref()); - } - - /** Optimized triangular matrix - matrix product */ - template - Eval > operator * (const MatrixBase& other) const - { - return this->lazyProduct(other).eval(); - } - - /** Optimized matrix - triangular matrix product */ - template - friend Eval, Transpose > > > - operator * (const MatrixBase& other, const Triangular& 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 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 +// +// 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_TRIANGULAR_ASSIGN_H +#define EIGEN_TRIANGULAR_ASSIGN_H + +template +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::run(dst, src); + dst.coeffRef(row, col) = src.coeff(row, col); + } +}; + +template +struct ei_triangular_assign_unroller +{ + 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 +struct ei_triangular_assign_unroller +{ + static void run(Derived1 &, const Derived2 &) {} +}; + +template +struct ei_triangular_assign_unroller +{ + static void run(Derived1 &, const Derived2 &) {} +}; + + +template +struct ei_assignment_impl +{ + 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 + ::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