aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/LU1
-rw-r--r--Eigen/src/Core/Block.h16
-rw-r--r--Eigen/src/Core/MatrixBase.h1
-rw-r--r--Eigen/src/LU/Determinant.h78
-rw-r--r--Eigen/src/LU/Inverse.h84
-rw-r--r--bench/benchmark.cpp2
6 files changed, 150 insertions, 32 deletions
diff --git a/Eigen/LU b/Eigen/LU
index 9d92196d5..c22eac2b5 100644
--- a/Eigen/LU
+++ b/Eigen/LU
@@ -5,6 +5,7 @@
namespace Eigen {
+#include "Eigen/src/LU/Determinant.h"
#include "Eigen/src/LU/Inverse.h"
} // namespace Eigen
diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h
index 0ff72075e..2c9f429a2 100644
--- a/Eigen/src/Core/Block.h
+++ b/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.");
}
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index d5912c9e6..23a82051f 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/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:
diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
new file mode 100644
index 000000000..16c33ec31
--- /dev/null
+++ b/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 \ No newline at end of file
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 4a0cd9f40..da74cc4a5 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/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++)
{
@@ -137,6 +142,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 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 if(_Size == 2)
{
- 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 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
+ else if(_Size == 3)
{
- _compute_not_unrolled(matrix, max);
+ 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.
diff --git a/bench/benchmark.cpp b/bench/benchmark.cpp
index 418bf72f2..39ed317aa 100644
--- a/bench/benchmark.cpp
+++ b/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;
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-27 04:36:57 +0000