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 MatrixBase { case TopLeft: return Block(derived(), 0, 0, cRows, cCols); - break; case TopRight: return Block(derived(), 0, cols() - cCols, cRows, cCols); - break; case BottomLeft: return Block(derived(), rows() - cRows, 0, cRows, cCols); - break; case BottomRight: return Block(derived(), rows() - cRows, cols() - cCols, cRows, cCols); - break; default: ei_assert(false && "Bad corner type."); } @@ -416,16 +412,12 @@ const Block MatrixBase { case TopLeft: return Block(derived(), 0, 0, cRows, cCols); - break; case TopRight: return Block(derived(), 0, cols() - cCols, cRows, cCols); - break; case BottomLeft: return Block(derived(), rows() - cRows, 0, cRows, cCols); - break; case BottomRight: return Block(derived(), rows() - cRows, cols() - cCols, cRows, cCols); - break; default: ei_assert(false && "Bad corner type."); } @@ -452,16 +444,12 @@ Block MatrixBase { case TopLeft: return Block(derived(), 0, 0); - break; case TopRight: return Block(derived(), 0, cols() - CCols); - break; case BottomLeft: return Block(derived(), rows() - CRows, 0); - break; case BottomRight: return Block(derived(), rows() - CRows, cols() - CCols); - break; default: ei_assert(false && "Bad corner type."); } @@ -477,16 +465,12 @@ const Block MatrixBase { case TopLeft: return Block(derived(), 0, 0); - break; case TopRight: return Block(derived(), 0, cols() - CCols); - break; case BottomLeft: return Block(derived(), rows() - CRows, 0); - break; case BottomRight: return Block(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 class MatrixBase //@{ const Inverse inverse() const; const Inverse 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 +// +// 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 . + +#ifndef EIGEN_DETERMINANT_H +#define EIGEN_DETERMINANT_H + +template +const typename Derived::Scalar ei_bruteforce_det3_helper +(const MatrixBase& 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 +const typename Derived::Scalar ei_bruteforce_det4_helper +(const MatrixBase& 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 +const typename Derived::Scalar ei_bruteforce_det(const MatrixBase& 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 ei_traits::Scalar MatrixBase::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 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 struct _unroll_first_loop; template struct _unroll_second_loop; template struct _unroll_third_loop; @@ -99,8 +101,11 @@ template class Inverse : ei_no_ass template void Inverse -::_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(0); const int size = matrix.rows(); for(int k = 0; k < size-1; k++) { @@ -136,6 +141,21 @@ void Inverse } } +template +void Inverse +::_compute_unrolled(const ExpressionType& xpr) +{ + MatrixType matrix(xpr); + const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff() + : static_cast(0); + const int size = MatrixType::RowsAtCompileTime; + _unroll_first_loop::run(*this, matrix, max); + if(CheckExistence && !m_exists) return; + _unroll_second_loop::run(*this, matrix, max); + if(CheckExistence && !m_exists) return; + _unroll_third_loop::run(*this, matrix, max); +} + template template struct Inverse::_unroll_first_loop @@ -246,23 +266,57 @@ struct Inverse::_unroll_third_loop void Inverse::_compute(const ExpressionType& xpr) { - MatrixType matrix(xpr); - const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff() - : static_cast(0); - - if(MatrixType::RowsAtCompileTime <= 4) + if(_Size == 1) { - const int size = MatrixType::RowsAtCompileTime; - _unroll_first_loop::run(*this, matrix, max); - if(CheckExistence && !m_exists) return; - _unroll_second_loop::run(*this, matrix, max); - if(CheckExistence && !m_exists) return; - _unroll_third_loop::run(*this, matrix, max); + const Scalar x = xpr.coeff(0,0); + if(CheckExistence && x == static_cast(0)) + m_exists = false; + else + m_inverse.coeffRef(0,0) = static_cast(1) / x; } - else + else if(_Size == 2) { - _compute_not_unrolled(matrix, max); + const typename ei_nested::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(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::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(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::ones() + 0.00005 * (m + m*m); + m = I + 0.00005 * (m + m*m); } asm("#end"); cout << m << endl;