add dot product, extend meta unrolling everywhere

src/Core.h

@@ -12,3 +12,4 @@
 #include "Core/Transpose.h"
 #include "Core/Conjugate.h"
 #include "Core/Trace.h"
+#include "Core/Dot.h"

src/Core/CopyHelper.h

@@ -27,37 +27,44 @@
 #ifndef EI_COPYHELPER_H
 #define EI_COPYHELPER_H
 
-template<int UnrollCount, int Rows> class EiCopyHelperUnroller
+template<int UnrollCount, int Rows> struct EiCopyHelperUnroller
 {
-    static const int col = (UnrollCount-1) / Rows;
+  static const int col = (UnrollCount-1) / Rows;
-    static const int row = (UnrollCount-1) % Rows;
+  static const int row = (UnrollCount-1) % Rows;
 
-  public:
+  template <typename Derived1, typename Derived2>
-    template <typename Derived1, typename Derived2>
+  static void run(Derived1 &dst, const Derived2 &src)
-    static void run(Derived1 &dst, const Derived2 &src)
+  {
-    {
+    EiCopyHelperUnroller<UnrollCount-1, Rows>::run(dst, src);
-      EiCopyHelperUnroller<UnrollCount-1, Rows>::run(dst, src);
+    dst.write(row, col) = src.read(row, col);
-      dst.write(row, col) = src.read(row, col);
+  }
-    }
 };
 
-template<int Rows> class EiCopyHelperUnroller<0, Rows>
+template<int Rows> struct EiCopyHelperUnroller<0, Rows>
 {
-  public:
+  template <typename Derived1, typename Derived2>
-    template <typename Derived1, typename Derived2>
+  static void run(Derived1 &dst, const Derived2 &src)
-    static void run(Derived1 &dst, const Derived2 &src)
+  {
-    {
+    dst.write(0, 0) = src.read(0, 0);
-      EI_UNUSED(dst);
+  }
-      EI_UNUSED(src);
+};
-    }
+
+template<int Rows> struct EiCopyHelperUnroller<EiDynamic, Rows>
+{
+  template <typename Derived1, typename Derived2>
+  static void run(Derived1 &dst, const Derived2 &src)
+  {
+    EI_UNUSED(dst);
+    EI_UNUSED(src);
+  }
 };
 
 template<typename Scalar, typename Derived>
 template<typename OtherDerived>
 void EiObject<Scalar, Derived>::_copy_helper(const EiObject<Scalar, OtherDerived>& other)
 {
-  if(UnrollCount > 0 && UnrollCount <= EI_LOOP_UNROLLING_LIMIT)
+  if(SizeAtCompileTime != EiDynamic && SizeAtCompileTime <= EI_LOOP_UNROLLING_LIMIT)
-    EiCopyHelperUnroller<UnrollCount, RowsAtCompileTime>::run(*this, other);
+    EiCopyHelperUnroller<SizeAtCompileTime, RowsAtCompileTime>::run(*this, other);
   else
     for(int i = 0; i < rows(); i++)
       for(int j = 0; j < cols(); j++)

src/Core/Dot.h

@@ -0,0 +1,83 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2007 Benoit Jacob <jacob@math.jussieu.fr>
+//
+// Eigen is free software; 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 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 General Public License for more
+// details.
+//
+// You should have received a copy of the GNU General Public License along
+// with Eigen; if not, write to the Free Software Foundation, Inc., 51
+// Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
+//
+// As a special exception, if other files instantiate templates or use macros
+// or functions from this file, or you compile this file and link it
+// with other works to produce a work based on this file, this file does not
+// by itself cause the resulting work to be covered by the GNU General Public
+// License. This exception does not invalidate any other reasons why a work
+// based on this file might be covered by the GNU General Public License.
+
+#ifndef EI_DOT_H
+#define EI_DOT_H
+
+template<int Index, int Size, typename Derived1, typename Derived2>
+struct EiDotUnroller
+{
+  static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot)
+  {
+    const int i = Index - 1;
+    if(i == Size - 1)
+      dot = v1[i] * EiConj(v2[i]);
+    else
+      dot += v1[i] * EiConj(v2[i]);
+    EiDotUnroller<Index-1, Size, Derived1, Derived2>::run(v1, v2, dot);
+  }
+};
+
+template<int Size, typename Derived1, typename Derived2>
+struct EiDotUnroller<0, Size, Derived1, Derived2>
+{
+  static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot)
+  {
+    EI_UNUSED(v1);
+    EI_UNUSED(v2);
+    EI_UNUSED(dot);
+  }
+};
+
+template<int Size, typename Derived1, typename Derived2>
+struct EiDotUnroller<EiDynamic, Size, Derived1, Derived2>
+{
+  static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot)
+  {
+    EI_UNUSED(v1);
+    EI_UNUSED(v2);
+    EI_UNUSED(dot);
+  }
+};
+
+template<typename Scalar, typename Derived>
+template<typename OtherDerived>
+Scalar EiObject<Scalar, Derived>::dot(const OtherDerived& other) const
+{
+  assert(IsVector && OtherDerived::IsVector && size() == other.size());
+  Scalar res;
+  if(SizeAtCompileTime != EiDynamic && SizeAtCompileTime <= 16)
+    EiDotUnroller<SizeAtCompileTime, SizeAtCompileTime, Derived, OtherDerived>
+      ::run(*static_cast<const Derived*>(this), other, res);
+  else
+  {
+    res = (*this)[0] * EiConj(other[0]);
+    for(int i = 1; i < size(); i++)
+      res += (*this)[i]* EiConj(other[i]);
+  }
+  return res;
+}
+
+#endif // EI_DOT_H

src/Core/MatrixOps.h

@@ -105,6 +105,49 @@ template<typename Lhs, typename Rhs> class EiDifference
     const RhsRef m_rhs;
 };
 
+template<int Index, int Size, typename Lhs, typename Rhs>
+struct EiMatrixProductUnroller
+{
+  static void run(int row, int col, const Lhs& lhs, const Rhs& rhs,
+                  typename Lhs::Scalar &res)
+  {
+    const int i = Index - 1;
+    if(i == Size - 1)
+      res = lhs.read(row, i) * rhs.read(i, col);
+    else
+      res += lhs.read(row, i) * rhs.read(i, col);
+    EiMatrixProductUnroller<Index-1, Size, Lhs, Rhs>::run(row, col, lhs, rhs, res);
+  }
+};
+
+template<int Size, typename Lhs, typename Rhs>
+struct EiMatrixProductUnroller<0, Size, Lhs, Rhs>
+{
+  static void run(int row, int col, const Lhs& lhs, const Rhs& rhs,
+                  typename Lhs::Scalar &res)
+  {
+    EI_UNUSED(row);
+    EI_UNUSED(col);
+    EI_UNUSED(lhs);
+    EI_UNUSED(rhs);
+    EI_UNUSED(res);
+  }
+};
+
+template<int Size, typename Lhs, typename Rhs>
+struct EiMatrixProductUnroller<EiDynamic, Size, Lhs, Rhs>
+{
+  static void run(int row, int col, const Lhs& lhs, const Rhs& rhs,
+                  typename Lhs::Scalar &res)
+  {
+    EI_UNUSED(row);
+    EI_UNUSED(col);
+    EI_UNUSED(lhs);
+    EI_UNUSED(rhs);
+    EI_UNUSED(res);
+  }
+};
+
 template<typename Lhs, typename Rhs> class EiMatrixProduct
   : public EiObject<typename Lhs::Scalar, EiMatrixProduct<Lhs, Rhs> >
 {
@@ -136,17 +179,17 @@ template<typename Lhs, typename Rhs> class EiMatrixProduct
     
     Scalar _read(int row, int col) const
     {
-      if(Lhs::ColsAtCompileTime == 3)
+      Scalar res;
-      {
+      if(Lhs::ColsAtCompileTime != EiDynamic && Lhs::ColsAtCompileTime <= 16)
-        return m_lhs(row,0) * m_rhs(0,col) + m_lhs(row,1) * m_rhs(1,col) + m_lhs(row,2) * m_rhs(2,col);
+        EiMatrixProductUnroller<Lhs::ColsAtCompileTime, Lhs::ColsAtCompileTime, LhsRef, RhsRef>
-      }
+          ::run(row, col, m_lhs, m_rhs, res);
       else
       {
-        Scalar x = static_cast<Scalar>(0);
+        res = m_lhs(row, 0) * m_rhs(0, col);
-        for(int i = 0; i < m_lhs.cols(); i++)
+        for(int i = 1; i < m_lhs.cols(); i++)
-          x += m_lhs.read(row, i) * m_rhs.read(i, col);
+          res += m_lhs(row, i) * m_rhs(i, col);
-        return x;
       }
+      return res;
     }
     
   protected:

src/Core/Object.h

@@ -29,11 +29,11 @@
 template<typename Scalar, typename Derived> class EiObject
 {
     static const int RowsAtCompileTime = Derived::RowsAtCompileTime,
-                     ColsAtCompileTime = Derived::ColsAtCompileTime;
+                     ColsAtCompileTime = Derived::ColsAtCompileTime,
-    static const bool HasDynamicSize = RowsAtCompileTime != EiDynamic
+                     SizeAtCompileTime
-                                    && ColsAtCompileTime != EiDynamic;
+      = RowsAtCompileTime == EiDynamic || ColsAtCompileTime == EiDynamic
-    static const int UnrollCount = HasDynamicSize ?
+      ? EiDynamic : RowsAtCompileTime * ColsAtCompileTime;
-                                   RowsAtCompileTime * ColsAtCompileTime : 0;
+    static const bool IsVector = RowsAtCompileTime == 1 || ColsAtCompileTime == 1;
     
     template<typename OtherDerived>
     void _copy_helper(const EiObject<Scalar, OtherDerived>& other);
@@ -88,6 +88,9 @@ template<typename Scalar, typename Derived> class EiObject
     EiTranspose<EiConjugate<Derived> > adjoint() const { return conjugate().transpose(); }
     Scalar trace() const;
     
+    template<typename OtherDerived>
+    Scalar dot(const OtherDerived& other) const;
+    
     template<typename OtherDerived>
     EiMatrixProduct<Derived, OtherDerived>
     lazyMul(const EiObject<Scalar, OtherDerived>& other) const EI_ALWAYS_INLINE;
@@ -121,14 +124,14 @@ template<typename Scalar, typename Derived> class EiObject
     
     Scalar operator[](int index) const
     {
-      assert(RowsAtCompileTime == 1 || ColsAtCompileTime == 1);
+      assert(IsVector);
       if(RowsAtCompileTime == 1) return read(0, index);
       else return read(index, 0);
     }
     
     Scalar& operator[](int index)
     {
-      assert(RowsAtCompileTime == 1 || ColsAtCompileTime == 1);
+      assert(IsVector);
       if(RowsAtCompileTime == 1) return write(0, index);
       else return write(index, 0);
     }

src/Core/Trace.h

@@ -1,7 +1,6 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra. Eigen itself is part of the KDE project.
 //
-// Copyright (C) 2007 Michael Olbrich <michael.olbrich@gmx.net>
 // Copyright (C) 2006-2007 Benoit Jacob <jacob@math.jussieu.fr>
 //
 // Eigen is free software; you can redistribute it and/or modify it under the
@@ -27,25 +26,31 @@
 #ifndef EI_TRACE_H
 #define EI_TRACE_H
 
-template<int CurrentRow, int Rows, typename Derived> struct EiTraceUnroller
+template<int Index, int Rows, typename Derived> struct EiTraceUnroller
 {
-  typedef typename Derived::Scalar Scalar;
+  static void run(const Derived &mat, typename Derived::Scalar &trace)
-
-  static void run(const Derived &mat, Scalar &trace)
   {
-    if(CurrentRow == Rows - 1)
+    const int i = Index - 1;
-      trace = mat(CurrentRow, CurrentRow);
+    if(i == Rows - 1)
+      trace = mat(i, i);
     else
-      trace += mat(CurrentRow, CurrentRow);
+      trace += mat(i, i);
-    EiTraceUnroller<CurrentRow-1, Rows, Derived>::run(mat, trace);
+    EiTraceUnroller<Index-1, Rows, Derived>::run(mat, trace);
   }
 };
 
-template<int Rows, typename Derived> struct EiTraceUnroller<-1, Rows, Derived>
+template<int Rows, typename Derived> struct EiTraceUnroller<0, Rows, Derived>
 {
-  typedef typename Derived::Scalar Scalar;
+  static void run(const Derived &mat, typename Derived::Scalar &trace)
+  {
+    EI_UNUSED(mat);
+    EI_UNUSED(trace);
+  }
+};
 
-  static void run(const Derived &mat, Scalar &trace)
+template<int Rows, typename Derived> struct EiTraceUnroller<EiDynamic, Rows, Derived>
+{
+  static void run(const Derived &mat, typename Derived::Scalar &trace)
   {
     EI_UNUSED(mat);
     EI_UNUSED(trace);
@@ -58,7 +63,7 @@ Scalar EiObject<Scalar, Derived>::trace() const
   assert(rows() == cols());
   Scalar res;
   if(RowsAtCompileTime != EiDynamic && RowsAtCompileTime <= 16)
-    EiTraceUnroller<RowsAtCompileTime - 1, RowsAtCompileTime, Derived>
+    EiTraceUnroller<RowsAtCompileTime, RowsAtCompileTime, Derived>
       ::run(*static_cast<const Derived*>(this), res);
   else
   {