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

@@ -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
 

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();
 }
 

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);
           }
         }
       }

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