- optimized determinant calculations for small matrices (size <= 4)

(only 30 muls for size 4)
- rework the matrix inversion: now using cofactor technique for size<=3,
  so the ugly unrolling is only used for size 4 anymore, and even there
  I'm looking to get rid of it.

Eigen/LU

@@ -5,6 +5,7 @@
 
 namespace Eigen {
 
+#include "Eigen/src/LU/Determinant.h"
 #include "Eigen/src/LU/Inverse.h"
 
 } // namespace Eigen

Eigen/src/Core/Block.h

@@ -392,16 +392,12 @@ Block<Derived> MatrixBase<Derived>
   {
     case TopLeft:
       return Block<Derived>(derived(), 0, 0, cRows, cCols);
-      break;
     case TopRight:
       return Block<Derived>(derived(), 0, cols() - cCols, cRows, cCols);
-      break;
     case BottomLeft:
       return Block<Derived>(derived(), rows() - cRows, 0, cRows, cCols);
-      break;
     case BottomRight:
       return Block<Derived>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
-      break;
     default:
       ei_assert(false && "Bad corner type.");
   }
@@ -416,16 +412,12 @@ const Block<Derived> MatrixBase<Derived>
   {
     case TopLeft:
       return Block<Derived>(derived(), 0, 0, cRows, cCols);
-      break;
     case TopRight:
       return Block<Derived>(derived(), 0, cols() - cCols, cRows, cCols);
-      break;
     case BottomLeft:
       return Block<Derived>(derived(), rows() - cRows, 0, cRows, cCols);
-      break;
     case BottomRight:
       return Block<Derived>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
-      break;
     default:
       ei_assert(false && "Bad corner type.");
   }
@@ -452,16 +444,12 @@ Block<Derived, CRows, CCols> MatrixBase<Derived>
   {
     case TopLeft:
       return Block<Derived, CRows, CCols>(derived(), 0, 0);
-      break;
     case TopRight:
       return Block<Derived, CRows, CCols>(derived(), 0, cols() - CCols);
-      break;
     case BottomLeft:
       return Block<Derived, CRows, CCols>(derived(), rows() - CRows, 0);
-      break;
     case BottomRight:
       return Block<Derived, CRows, CCols>(derived(), rows() - CRows, cols() - CCols);
-      break;
     default:
       ei_assert(false && "Bad corner type.");
   }
@@ -477,16 +465,12 @@ const Block<Derived, CRows, CCols> MatrixBase<Derived>
   {
     case TopLeft:
       return Block<Derived, CRows, CCols>(derived(), 0, 0);
-      break;
     case TopRight:
       return Block<Derived, CRows, CCols>(derived(), 0, cols() - CCols);
-      break;
     case BottomLeft:
       return Block<Derived, CRows, CCols>(derived(), rows() - CRows, 0);
-      break;
     case BottomRight:
       return Block<Derived, CRows, CCols>(derived(), rows() - CRows, cols() - CCols);
-      break;
     default:
       ei_assert(false && "Bad corner type.");
   }

Eigen/src/Core/MatrixBase.h

@@ -513,6 +513,7 @@ template<typename Derived> class MatrixBase
     //@{
     const Inverse<Derived, true> inverse() const;
     const Inverse<Derived, false> quickInverse() const;
+    Scalar determinant() const;
     //@}
 
   private:

Eigen/src/LU/Determinant.h

@@ -0,0 +1,78 @@
+// 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 Benoit Jacob <jacob@math.jussieu.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_DETERMINANT_H
+#define EIGEN_DETERMINANT_H
+
+template<typename Derived>
+const typename Derived::Scalar ei_bruteforce_det3_helper
+(const MatrixBase<Derived>& matrix, int a, int b, int c)
+{
+  return matrix.coeff(0,a)
+         * (matrix.coeff(1,b) * matrix.coeff(2,c) - matrix.coeff(1,c) * matrix.coeff(2,b));
+}
+
+template<typename Derived>
+const typename Derived::Scalar ei_bruteforce_det4_helper
+(const MatrixBase<Derived>& matrix, int j, int k, int m, int n)
+{
+  return (matrix.coeff(j,0) * matrix.coeff(k,1) - matrix.coeff(k,0) * matrix.coeff(j,1))
+       * (matrix.coeff(m,2) * matrix.coeff(n,3) - matrix.coeff(n,2) * matrix.coeff(m,3));
+}
+
+template<typename Derived>
+const typename Derived::Scalar ei_bruteforce_det(const MatrixBase<Derived>& m)
+{
+  switch(Derived::RowsAtCompileTime)
+  {
+    case 1:
+      return m.coeff(0,0);
+    case 2:
+      return m.coeff(0,0) * m.coeff(1,1) - m.coeff(1,0) * m.coeff(0,1);
+    case 3:
+      return ei_bruteforce_det3_helper(m,0,1,2)
+           - ei_bruteforce_det3_helper(m,1,0,2)
+           + ei_bruteforce_det3_helper(m,2,0,1);
+    case 4:
+      // trick by Martin Costabel to compute 4x4 det with only 30 muls
+      return ei_bruteforce_det4_helper(m,0,1,2,3)
+           + ei_bruteforce_det4_helper(m,0,2,1,3)
+           + ei_bruteforce_det4_helper(m,0,3,1,2)
+           + ei_bruteforce_det4_helper(m,1,2,0,3)
+           + ei_bruteforce_det4_helper(m,1,3,0,2)
+           + ei_bruteforce_det4_helper(m,2,3,0,1);
+    default:
+      assert(false);
+  }
+}
+
+template<typename Derived>
+typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
+{
+  assert(rows() == cols());
+  if(rows() <= 4) return ei_bruteforce_det(derived());
+  else assert(false); // unimplemented for now
+}
+
+#endif // EIGEN_DETERMINANT_H

Eigen/src/LU/Inverse.h

@@ -86,8 +86,10 @@ template<typename ExpressionType, bool CheckExistence> class Inverse : ei_no_ass
       return m_inverse.coeff(row, col);
     }
 
+    enum { _Size = MatrixType::RowsAtCompileTime };
     void _compute(const ExpressionType& xpr);
-    void _compute_not_unrolled(MatrixType& matrix, const RealScalar& max);
+    void _compute_in_general_case(const ExpressionType& xpr);
+    void _compute_unrolled(const ExpressionType& xpr);
     template<int Size, int Step, bool Finished = Size==Dynamic> struct _unroll_first_loop;
     template<int Size, int Step, bool Finished = Size==Dynamic> struct _unroll_second_loop;
     template<int Size, int Step, bool Finished = Size==Dynamic> struct _unroll_third_loop;
@@ -99,8 +101,11 @@ template<typename ExpressionType, bool CheckExistence> class Inverse : ei_no_ass
 
 template<typename ExpressionType, bool CheckExistence>
 void Inverse<ExpressionType, CheckExistence>
-::_compute_not_unrolled(MatrixType& matrix, const RealScalar& max)
+::_compute_in_general_case(const ExpressionType& xpr)
 {
+  MatrixType matrix(xpr);
+  const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff()
+                                        : static_cast<RealScalar>(0);
   const int size = matrix.rows();
   for(int k = 0; k < size-1; k++)
   {
@@ -136,6 +141,21 @@ void Inverse<ExpressionType, CheckExistence>
   }
 }
 
+template<typename ExpressionType, bool CheckExistence>
+void Inverse<ExpressionType, CheckExistence>
+::_compute_unrolled(const ExpressionType& xpr)
+{
+  MatrixType matrix(xpr);
+  const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff()
+                                        : static_cast<RealScalar>(0);
+  const int size = MatrixType::RowsAtCompileTime;
+  _unroll_first_loop<size, 0>::run(*this, matrix, max);
+  if(CheckExistence && !m_exists) return;
+  _unroll_second_loop<size, 0>::run(*this, matrix, max);
+  if(CheckExistence && !m_exists) return;
+  _unroll_third_loop<size, 1>::run(*this, matrix, max);
+}
+
 template<typename ExpressionType, bool CheckExistence>
 template<int Size, int Step, bool Finished>
 struct Inverse<ExpressionType, CheckExistence>::_unroll_first_loop
@@ -246,23 +266,57 @@ struct Inverse<ExpressionType, CheckExistence>::_unroll_third_loop<Step, Size, t
 template<typename ExpressionType, bool CheckExistence>
 void Inverse<ExpressionType, CheckExistence>::_compute(const ExpressionType& xpr)
 {
-  MatrixType matrix(xpr);
-  const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff()
-                                        : static_cast<RealScalar>(0);
-
-  if(MatrixType::RowsAtCompileTime <= 4)
+  if(_Size == 1)
   {
-    const int size = MatrixType::RowsAtCompileTime;
-    _unroll_first_loop<size, 0>::run(*this, matrix, max);
-    if(CheckExistence && !m_exists) return;
-    _unroll_second_loop<size, 0>::run(*this, matrix, max);
-    if(CheckExistence && !m_exists) return;
-    _unroll_third_loop<size, 1>::run(*this, matrix, max);
+    const Scalar x = xpr.coeff(0,0);
+    if(CheckExistence && x == static_cast<Scalar>(0))
+      m_exists = false;
+    else
+      m_inverse.coeffRef(0,0) = static_cast<Scalar>(1) / x;
   }
-  else
+  else if(_Size == 2)
   {
-    _compute_not_unrolled(matrix, max);
+    const typename ei_nested<ExpressionType, 1+CheckExistence>::type matrix(xpr);
+    const Scalar det = matrix.determinant();
+    if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwiseAbs().maxCoeff()))
+      m_exists = false;
+    else
+    {
+      const Scalar invdet = static_cast<Scalar>(1) / det;
+      m_inverse.coeffRef(0,0) = matrix.coeff(1,1) * invdet;
+      m_inverse.coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
+      m_inverse.coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
+      m_inverse.coeffRef(1,1) = matrix.coeff(0,0) * invdet;
+    }
   }
+  else if(_Size == 3)
+  {
+    const typename ei_nested<ExpressionType, 2+CheckExistence>::type matrix(xpr);
+    const Scalar det_minor00 = matrix.minor(0,0).determinant();
+    const Scalar det_minor10 = matrix.minor(1,0).determinant();
+    const Scalar det_minor20 = matrix.minor(2,0).determinant();
+    const Scalar det = det_minor00 * matrix.coeff(0,0)
+                      - det_minor10 * matrix.coeff(1,0)
+                      + det_minor20 * matrix.coeff(2,0);
+    if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwiseAbs().maxCoeff()))
+      m_exists = false;
+    else
+    {
+      const Scalar invdet = static_cast<Scalar>(1) / det;
+      m_inverse.coeffRef(0, 0) = det_minor00 * invdet;
+      m_inverse.coeffRef(0, 1) = -det_minor10 * invdet;
+      m_inverse.coeffRef(0, 2) = det_minor20 * invdet;
+      m_inverse.coeffRef(1, 0) = -matrix.minor(0,1).determinant() * invdet;
+      m_inverse.coeffRef(1, 1) = matrix.minor(1,1).determinant() * invdet;
+      m_inverse.coeffRef(1, 2) = -matrix.minor(2,1).determinant() * invdet;
+      m_inverse.coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
+      m_inverse.coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
+      m_inverse.coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
+    }
+  }
+  else if(_Size == 4)
+    _compute_unrolled(xpr);
+  else _compute_in_general_case(xpr);
 }
 
 /** \return the matrix inverse of \c *this, if it exists.

bench/benchmark.cpp

@@ -28,7 +28,7 @@ int main(int argc, char *argv[])
     asm("#begin");
     for(int a = 0; a < REPEAT; a++)
     {
-        m = Matrix<SCALAR,MATSIZE,MATSIZE>::ones() + 0.00005 * (m + m*m);
+        m = I + 0.00005 * (m + m*m);
     }
     asm("#end");
     cout << m << endl;