diff options
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/LU | 1 | ||||
-rw-r--r-- | Eigen/src/Core/Block.h | 16 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 1 | ||||
-rw-r--r-- | Eigen/src/LU/Determinant.h | 78 | ||||
-rw-r--r-- | Eigen/src/LU/Inverse.h | 84 |
5 files changed, 149 insertions, 31 deletions
@@ -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. |