reorganize meta loop unrolling, add Trace (meta-unrolled), fix compilation issues in

the conjugation/adjunction stuff

doc/tutorial.cpp

@@ -30,6 +30,5 @@ int main(int, char **)
   cout << "Column 1 of m2 is:" << endl << m2.col(1) << endl;
   cout << "The transpose of m2 is:" << endl << m2.transpose() << endl;
   cout << "The matrix m2 with row 0 and column 1 removed is:" << endl << m2.minor(0,1) << endl;
-  
   return 0;
 }

src/Core.h

@@ -1,6 +1,7 @@
 #include "Core/Util.h"
 #include "Core/Numeric.h"
 #include "Core/Object.h"
+#include "Core/CopyHelper.h"
 #include "Core/MatrixRef.h"
 #include "Core/MatrixStorage.h"
 #include "Core/Matrix.h"
@@ -10,3 +11,4 @@
 #include "Core/Minor.h"
 #include "Core/Transpose.h"
 #include "Core/Conjugate.h"
+#include "Core/Trace.h"

src/Core/Conjugate.h

@@ -31,7 +31,7 @@ template<typename MatrixType> class EiConjugate
 {
   public:
     typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::Ref MatRef;
+    typedef typename MatrixType::ConstRef MatRef;
     friend class EiObject<Scalar, EiConjugate<MatrixType> >;
     
     static const int RowsAtCompileTime = MatrixType::RowsAtCompileTime,
@@ -61,9 +61,9 @@ template<typename MatrixType> class EiConjugate
 
 template<typename Scalar, typename Derived>
 EiConjugate<Derived>
-EiObject<Scalar, Derived>::conjugate()
+EiObject<Scalar, Derived>::conjugate() const
 {
-  return EiConjugate<Derived>(static_cast<Derived*>(this)->ref());
+  return EiConjugate<Derived>(static_cast<const Derived*>(this)->constRef());
 }
 
 #endif // EI_CONJUGATE_H

src/Core/CopyHelper.h

@@ -2,6 +2,7 @@
 // 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
 // terms of the GNU General Public License as published by the Free Software
@@ -23,32 +24,44 @@
 // 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_LOOP_H
-#define EI_LOOP_H
+#ifndef EI_COPYHELPER_H
+#define EI_COPYHELPER_H
 
-template<int UnrollCount, int Rows> class EiLoop
+template<int UnrollCount, int Rows> class EiCopyHelperUnroller
 {
     static const int col = (UnrollCount-1) / Rows;
     static const int row = (UnrollCount-1) % Rows;
 
   public:
     template <typename Derived1, typename Derived2>
-    static void copy(Derived1 &dst, const Derived2 &src)
+    static void run(Derived1 &dst, const Derived2 &src)
     {
-      EiLoop<UnrollCount-1, Rows>::copy(dst, src);
+      EiCopyHelperUnroller<UnrollCount-1, Rows>::run(dst, src);
       dst.write(row, col) = src.read(row, col);
     }
 };
 
-template<int Rows> class EiLoop<0, Rows>
+template<int Rows> class EiCopyHelperUnroller<0, Rows>
 {
   public:
     template <typename Derived1, typename Derived2>
-    static void copy(Derived1 &dst, const Derived2 &src)
+    static void run(Derived1 &dst, const Derived2 &src)
     {
       EI_UNUSED(dst);
       EI_UNUSED(src);
     }
 };
 
-#endif // EI_LOOP_H
+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)
+    EiCopyHelperUnroller<UnrollCount, RowsAtCompileTime>::run(*this, other);
+  else
+    for(int i = 0; i < rows(); i++)
+      for(int j = 0; j < cols(); j++)
+        write(i, j) = other.read(i, j);
+}
+
+#endif // EI_COPYHELPER_H

src/Core/Object.h

@@ -36,15 +36,7 @@ template<typename Scalar, typename Derived> class EiObject
                                    RowsAtCompileTime * ColsAtCompileTime : 0;
     
     template<typename OtherDerived>
-    void _copy_helper(const EiObject<Scalar, OtherDerived>& other)
-    {
-      if(UnrollCount > 0 && UnrollCount <= EI_LOOP_UNROLLING_LIMIT)
-        EiLoop<UnrollCount, RowsAtCompileTime>::copy(*this, other);
-      else
-        for(int i = 0; i < rows(); i++)
-          for(int j = 0; j < cols(); j++)
-            write(i, j) = other.read(i, j);
-    }
+    void _copy_helper(const EiObject<Scalar, OtherDerived>& other);
     
   public:
     typedef typename EiForwardDecl<Derived>::Ref Ref;
@@ -92,8 +84,9 @@ template<typename Scalar, typename Derived> class EiObject
     EiMinor<Derived> minor(int row, int col);
     EiBlock<Derived> block(int startRow, int endRow, int startCol, int endCol);
     EiTranspose<Derived> transpose();
-    EiConjugate<Derived> conjugate();
-    EiTranspose<EiConjugate<Derived> > adjoint() { return conjugate().transpose(); }
+    EiConjugate<Derived> conjugate() const;
+    EiTranspose<EiConjugate<Derived> > adjoint() const { return conjugate().transpose(); }
+    Scalar trace() const;
     
     template<typename OtherDerived>
     EiMatrixProduct<Derived, OtherDerived>

src/Core/Trace.h

@@ -0,0 +1,72 @@
+// 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
+// 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_TRACE_H
+#define EI_TRACE_H
+
+template<int CurrentRow, int Rows, typename Derived> struct EiTraceUnroller
+{
+  typedef typename Derived::Scalar Scalar;
+
+  static void run(const Derived &mat, Scalar &trace)
+  {
+    if(CurrentRow == Rows - 1)
+      trace = mat(CurrentRow, CurrentRow);
+    else
+      trace += mat(CurrentRow, CurrentRow);
+    EiTraceUnroller<CurrentRow-1, Rows, Derived>::run(mat, trace);
+  }
+};
+
+template<int Rows, typename Derived> struct EiTraceUnroller<-1, Rows, Derived>
+{
+  typedef typename Derived::Scalar Scalar;
+
+  static void run(const Derived &mat, Scalar &trace)
+  {
+    EI_UNUSED(mat);
+    EI_UNUSED(trace);
+  }
+};
+
+template<typename Scalar, typename Derived>
+Scalar EiObject<Scalar, Derived>::trace() const
+{
+  assert(rows() == cols());
+  Scalar res;
+  if(RowsAtCompileTime != EiDynamic && RowsAtCompileTime <= 16)
+    EiTraceUnroller<RowsAtCompileTime - 1, RowsAtCompileTime, Derived>
+      ::run(*static_cast<const Derived*>(this), res);
+  else
+  {
+    res = read(0, 0);
+    for(int i = 1; i < rows(); i++)
+      res += read(i, i);
+  }
+  return res;
+}
+
+#endif // EI_TRACE_H