diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2008-04-14 17:07:12 +0000 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2008-04-14 17:07:12 +0000 |
commit | 2a86f052a5d26c6917aee91004e220c0c2a0657e (patch) | |
tree | 910b459e63c41ac0d66a8c39332c977029615d8c /Eigen/src/LU/Determinant.h | |
parent | 9789c04467d17abc338981cd2aa6d8824e6705b4 (diff) |
- 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.
Diffstat (limited to 'Eigen/src/LU/Determinant.h')
-rw-r--r-- | Eigen/src/LU/Determinant.h | 78 |
1 files changed, 78 insertions, 0 deletions
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 |