mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-08-12 03:39:01 +08:00
remove Minor
adapt 3x3 and 4x4 (non-SSE) inverse paths
This commit is contained in:
Benoit Jacob
parent
abbe260905
commit
2362eadcdd
@ -77,7 +77,7 @@
|
||||
#define EIGEN_VECTORIZE_SSE2
|
||||
|
||||
// Detect sse3/ssse3/sse4:
|
||||
// gcc and icc defines __SSE3__, ..,
|
||||
// gcc and icc defines __SSE3__, ...
|
||||
// there is no way to know about this on msvc. You can define EIGEN_VECTORIZE_SSE* if you
|
||||
// want to force the use of those instructions with msvc.
|
||||
#ifdef __SSE3__
|
||||
@ -169,9 +169,6 @@
|
||||
// defined in bits/termios.h
|
||||
#undef B0
|
||||
|
||||
// defined in some GNU standard header
|
||||
#undef minor
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
inline static const char *SimdInstructionSetsInUse(void) {
|
||||
@ -262,7 +259,6 @@ using std::size_t;
|
||||
#include "src/Core/Map.h"
|
||||
#include "src/Core/Block.h"
|
||||
#include "src/Core/VectorBlock.h"
|
||||
#include "src/Core/Minor.h"
|
||||
#include "src/Core/Transpose.h"
|
||||
#include "src/Core/DiagonalMatrix.h"
|
||||
#include "src/Core/Diagonal.h"
|
||||
|
||||
@ -208,9 +208,6 @@ template<typename Derived> class MatrixBase
|
||||
const AdjointReturnType adjoint() const;
|
||||
void adjointInPlace();
|
||||
|
||||
Minor<Derived> minor(int row, int col);
|
||||
const Minor<Derived> minor(int row, int col) const;
|
||||
|
||||
Diagonal<Derived,0> diagonal();
|
||||
const Diagonal<Derived,0> diagonal() const;
|
||||
|
||||
|
||||
@ -1,124 +0,0 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
|
||||
//
|
||||
// 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_MINOR_H
|
||||
#define EIGEN_MINOR_H
|
||||
|
||||
/** \nonstableyet
|
||||
* \class Minor
|
||||
*
|
||||
* \brief Expression of a minor
|
||||
*
|
||||
* \param MatrixType the type of the object in which we are taking a minor
|
||||
*
|
||||
* This class represents an expression of a minor. It is the return
|
||||
* type of MatrixBase::minor() and most of the time this is the only way it
|
||||
* is used.
|
||||
*
|
||||
* \sa MatrixBase::minor()
|
||||
*/
|
||||
template<typename MatrixType>
|
||||
struct ei_traits<Minor<MatrixType> >
|
||||
: ei_traits<MatrixType>
|
||||
{
|
||||
typedef typename ei_nested<MatrixType>::type MatrixTypeNested;
|
||||
typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
|
||||
enum {
|
||||
RowsAtCompileTime = (MatrixType::RowsAtCompileTime != Dynamic) ?
|
||||
int(MatrixType::RowsAtCompileTime) - 1 : Dynamic,
|
||||
ColsAtCompileTime = (MatrixType::ColsAtCompileTime != Dynamic) ?
|
||||
int(MatrixType::ColsAtCompileTime) - 1 : Dynamic,
|
||||
MaxRowsAtCompileTime = (MatrixType::MaxRowsAtCompileTime != Dynamic) ?
|
||||
int(MatrixType::MaxRowsAtCompileTime) - 1 : Dynamic,
|
||||
MaxColsAtCompileTime = (MatrixType::MaxColsAtCompileTime != Dynamic) ?
|
||||
int(MatrixType::MaxColsAtCompileTime) - 1 : Dynamic,
|
||||
Flags = _MatrixTypeNested::Flags & HereditaryBits,
|
||||
CoeffReadCost = _MatrixTypeNested::CoeffReadCost // minor is used typically on tiny matrices,
|
||||
// where loops are unrolled and the 'if' evaluates at compile time
|
||||
};
|
||||
};
|
||||
|
||||
template<typename MatrixType> class Minor
|
||||
: public MatrixBase<Minor<MatrixType> >
|
||||
{
|
||||
public:
|
||||
|
||||
typedef MatrixBase<Minor> Base;
|
||||
EIGEN_DENSE_PUBLIC_INTERFACE(Minor)
|
||||
|
||||
inline Minor(const MatrixType& matrix,
|
||||
int row, int col)
|
||||
: m_matrix(matrix), m_row(row), m_col(col)
|
||||
{
|
||||
ei_assert(row >= 0 && row < matrix.rows()
|
||||
&& col >= 0 && col < matrix.cols());
|
||||
}
|
||||
|
||||
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Minor)
|
||||
|
||||
inline int rows() const { return m_matrix.rows() - 1; }
|
||||
inline int cols() const { return m_matrix.cols() - 1; }
|
||||
|
||||
inline Scalar& coeffRef(int row, int col)
|
||||
{
|
||||
return m_matrix.const_cast_derived().coeffRef(row + (row >= m_row), col + (col >= m_col));
|
||||
}
|
||||
|
||||
inline const Scalar coeff(int row, int col) const
|
||||
{
|
||||
return m_matrix.coeff(row + (row >= m_row), col + (col >= m_col));
|
||||
}
|
||||
|
||||
protected:
|
||||
const typename MatrixType::Nested m_matrix;
|
||||
const int m_row, m_col;
|
||||
};
|
||||
|
||||
/** \nonstableyet
|
||||
* \return an expression of the (\a row, \a col)-minor of *this,
|
||||
* i.e. an expression constructed from *this by removing the specified
|
||||
* row and column.
|
||||
*
|
||||
* Example: \include MatrixBase_minor.cpp
|
||||
* Output: \verbinclude MatrixBase_minor.out
|
||||
*
|
||||
* \sa class Minor
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline Minor<Derived>
|
||||
MatrixBase<Derived>::minor(int row, int col)
|
||||
{
|
||||
return Minor<Derived>(derived(), row, col);
|
||||
}
|
||||
|
||||
/** \nonstableyet
|
||||
* This is the const version of minor(). */
|
||||
template<typename Derived>
|
||||
inline const Minor<Derived>
|
||||
MatrixBase<Derived>::minor(int row, int col) const
|
||||
{
|
||||
return Minor<Derived>(derived(), row, col);
|
||||
}
|
||||
|
||||
#endif // EIGEN_MINOR_H
|
||||
@ -48,7 +48,6 @@ template<typename ExpressionType, template <typename> class StorageBase > class
|
||||
template<typename ExpressionType> class NestByValue;
|
||||
template<typename ExpressionType> class ForceAlignedAccess;
|
||||
template<typename ExpressionType> class SwapWrapper;
|
||||
template<typename MatrixType> class Minor;
|
||||
|
||||
// MSVC has a big bug: when the expression ei_traits<MatrixType>::Flags&DirectAccessBit ? HasDirectAccess : NoDirectAccess
|
||||
// is used as default template parameter value here, it gets mis-evaluated as just ei_traits<MatrixType>::Flags
|
||||
|
||||
@ -122,6 +122,19 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 2>
|
||||
*** Size 3 implementation ***
|
||||
****************************/
|
||||
|
||||
template<typename MatrixType, int i, int j>
|
||||
inline typename MatrixType::Scalar ei_3x3_cofactor(const MatrixType& m)
|
||||
{
|
||||
enum {
|
||||
i1 = (i+1) % 3,
|
||||
i2 = (i+2) % 3,
|
||||
j1 = (j+1) % 3,
|
||||
j2 = (j+2) % 3
|
||||
};
|
||||
return m.coeff(i1, j1) * m.coeff(i2, j2)
|
||||
- m.coeff(i1, j2) * m.coeff(i2, j1);
|
||||
}
|
||||
|
||||
template<typename MatrixType, typename ResultType>
|
||||
inline void ei_compute_inverse_size3_helper(
|
||||
const MatrixType& matrix,
|
||||
@ -130,12 +143,12 @@ inline void ei_compute_inverse_size3_helper(
|
||||
ResultType& result)
|
||||
{
|
||||
result.row(0) = cofactors_col0 * invdet;
|
||||
result.coeffRef(1,0) = -matrix.minor(0,1).determinant() * invdet;
|
||||
result.coeffRef(1,1) = matrix.minor(1,1).determinant() * invdet;
|
||||
result.coeffRef(1,2) = -matrix.minor(2,1).determinant() * invdet;
|
||||
result.coeffRef(2,0) = matrix.minor(0,2).determinant() * invdet;
|
||||
result.coeffRef(2,1) = -matrix.minor(1,2).determinant() * invdet;
|
||||
result.coeffRef(2,2) = matrix.minor(2,2).determinant() * invdet;
|
||||
result.coeffRef(1,0) = ei_3x3_cofactor<MatrixType,0,1>(matrix) * invdet;
|
||||
result.coeffRef(1,1) = ei_3x3_cofactor<MatrixType,1,1>(matrix) * invdet;
|
||||
result.coeffRef(1,2) = ei_3x3_cofactor<MatrixType,2,1>(matrix) * invdet;
|
||||
result.coeffRef(2,0) = ei_3x3_cofactor<MatrixType,0,2>(matrix) * invdet;
|
||||
result.coeffRef(2,1) = ei_3x3_cofactor<MatrixType,1,2>(matrix) * invdet;
|
||||
result.coeffRef(2,2) = ei_3x3_cofactor<MatrixType,2,2>(matrix) * invdet;
|
||||
}
|
||||
|
||||
template<typename MatrixType, typename ResultType>
|
||||
@ -145,9 +158,9 @@ struct ei_compute_inverse<MatrixType, ResultType, 3>
|
||||
{
|
||||
typedef typename ResultType::Scalar Scalar;
|
||||
Matrix<Scalar,3,1> cofactors_col0;
|
||||
cofactors_col0.coeffRef(0) = matrix.minor(0,0).determinant();
|
||||
cofactors_col0.coeffRef(1) = -matrix.minor(1,0).determinant();
|
||||
cofactors_col0.coeffRef(2) = matrix.minor(2,0).determinant();
|
||||
cofactors_col0.coeffRef(0) = ei_3x3_cofactor<MatrixType,0,0>(matrix);
|
||||
cofactors_col0.coeffRef(1) = ei_3x3_cofactor<MatrixType,1,0>(matrix);
|
||||
cofactors_col0.coeffRef(2) = ei_3x3_cofactor<MatrixType,2,0>(matrix);
|
||||
const Scalar det = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
|
||||
const Scalar invdet = Scalar(1) / det;
|
||||
ei_compute_inverse_size3_helper(matrix, invdet, cofactors_col0, result);
|
||||
@ -167,9 +180,9 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
|
||||
{
|
||||
typedef typename ResultType::Scalar Scalar;
|
||||
Matrix<Scalar,3,1> cofactors_col0;
|
||||
cofactors_col0.coeffRef(0) = matrix.minor(0,0).determinant();
|
||||
cofactors_col0.coeffRef(1) = -matrix.minor(1,0).determinant();
|
||||
cofactors_col0.coeffRef(2) = matrix.minor(2,0).determinant();
|
||||
cofactors_col0.coeffRef(0) = ei_3x3_cofactor<MatrixType,0,0>(matrix);
|
||||
cofactors_col0.coeffRef(1) = ei_3x3_cofactor<MatrixType,1,0>(matrix);
|
||||
cofactors_col0.coeffRef(2) = ei_3x3_cofactor<MatrixType,2,0>(matrix);
|
||||
determinant = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
|
||||
invertible = ei_abs(determinant) > absDeterminantThreshold;
|
||||
if(!invertible) return;
|
||||
@ -182,27 +195,51 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
|
||||
*** Size 4 implementation ***
|
||||
****************************/
|
||||
|
||||
template<typename Derived>
|
||||
inline const typename Derived::Scalar ei_general_det3_helper
|
||||
(const MatrixBase<Derived>& matrix, int i1, int i2, int i3, int j1, int j2, int j3)
|
||||
{
|
||||
return matrix.coeff(i1,j1)
|
||||
* (matrix.coeff(i2,j2) * matrix.coeff(i3,j3) - matrix.coeff(i2,j3) * matrix.coeff(i3,j2));
|
||||
}
|
||||
|
||||
template<typename MatrixType, int i, int j>
|
||||
inline typename MatrixType::Scalar ei_4x4_cofactor(const MatrixType& matrix)
|
||||
{
|
||||
enum {
|
||||
i1 = (i+1) % 4,
|
||||
i2 = (i+2) % 4,
|
||||
i3 = (i+3) % 4,
|
||||
j1 = (j+1) % 4,
|
||||
j2 = (j+2) % 4,
|
||||
j3 = (j+3) % 4
|
||||
};
|
||||
return ei_general_det3_helper(matrix, i1, i2, i3, j1, j2, j3)
|
||||
+ ei_general_det3_helper(matrix, i2, i3, i1, j1, j2, j3)
|
||||
+ ei_general_det3_helper(matrix, i3, i1, i2, j1, j2, j3);
|
||||
}
|
||||
|
||||
template<int Arch, typename Scalar, typename MatrixType, typename ResultType>
|
||||
struct ei_compute_inverse_size4
|
||||
{
|
||||
static void run(const MatrixType& matrix, ResultType& result)
|
||||
{
|
||||
result.coeffRef(0,0) = matrix.minor(0,0).determinant();
|
||||
result.coeffRef(1,0) = -matrix.minor(0,1).determinant();
|
||||
result.coeffRef(2,0) = matrix.minor(0,2).determinant();
|
||||
result.coeffRef(3,0) = -matrix.minor(0,3).determinant();
|
||||
result.coeffRef(0,2) = matrix.minor(2,0).determinant();
|
||||
result.coeffRef(1,2) = -matrix.minor(2,1).determinant();
|
||||
result.coeffRef(2,2) = matrix.minor(2,2).determinant();
|
||||
result.coeffRef(3,2) = -matrix.minor(2,3).determinant();
|
||||
result.coeffRef(0,1) = -matrix.minor(1,0).determinant();
|
||||
result.coeffRef(1,1) = matrix.minor(1,1).determinant();
|
||||
result.coeffRef(2,1) = -matrix.minor(1,2).determinant();
|
||||
result.coeffRef(3,1) = matrix.minor(1,3).determinant();
|
||||
result.coeffRef(0,3) = -matrix.minor(3,0).determinant();
|
||||
result.coeffRef(1,3) = matrix.minor(3,1).determinant();
|
||||
result.coeffRef(2,3) = -matrix.minor(3,2).determinant();
|
||||
result.coeffRef(3,3) = matrix.minor(3,3).determinant();
|
||||
result.coeffRef(0,0) = ei_4x4_cofactor<MatrixType,0,0>(matrix);
|
||||
result.coeffRef(1,0) = -ei_4x4_cofactor<MatrixType,0,1>(matrix);
|
||||
result.coeffRef(2,0) = ei_4x4_cofactor<MatrixType,0,2>(matrix);
|
||||
result.coeffRef(3,0) = -ei_4x4_cofactor<MatrixType,0,3>(matrix);
|
||||
result.coeffRef(0,2) = ei_4x4_cofactor<MatrixType,2,0>(matrix);
|
||||
result.coeffRef(1,2) = -ei_4x4_cofactor<MatrixType,2,1>(matrix);
|
||||
result.coeffRef(2,2) = ei_4x4_cofactor<MatrixType,2,2>(matrix);
|
||||
result.coeffRef(3,2) = -ei_4x4_cofactor<MatrixType,2,3>(matrix);
|
||||
result.coeffRef(0,1) = -ei_4x4_cofactor<MatrixType,1,0>(matrix);
|
||||
result.coeffRef(1,1) = ei_4x4_cofactor<MatrixType,1,1>(matrix);
|
||||
result.coeffRef(2,1) = -ei_4x4_cofactor<MatrixType,1,2>(matrix);
|
||||
result.coeffRef(3,1) = ei_4x4_cofactor<MatrixType,1,3>(matrix);
|
||||
result.coeffRef(0,3) = -ei_4x4_cofactor<MatrixType,3,0>(matrix);
|
||||
result.coeffRef(1,3) = ei_4x4_cofactor<MatrixType,3,1>(matrix);
|
||||
result.coeffRef(2,3) = -ei_4x4_cofactor<MatrixType,3,2>(matrix);
|
||||
result.coeffRef(3,3) = ei_4x4_cofactor<MatrixType,3,3>(matrix);
|
||||
result /= (matrix.col(0).cwiseProduct(result.row(0).transpose())).sum();
|
||||
}
|
||||
};
|
||||
|
||||
@ -24,38 +24,8 @@
|
||||
|
||||
#include "main.h"
|
||||
|
||||
// check minor separately in order to avoid the possible creation of a zero-sized
|
||||
// array. Comes from a compilation error with gcc-3.4 or gcc-4 with -ansi -pedantic.
|
||||
// Another solution would be to declare the array like this: T m_data[Size==0?1:Size]; in ei_matrix_storage
|
||||
// but this is probably not bad to raise such an error at compile time...
|
||||
template<typename Scalar, int _Rows, int _Cols> struct CheckMinor
|
||||
{
|
||||
typedef Matrix<Scalar, _Rows, _Cols> MatrixType;
|
||||
CheckMinor(MatrixType& m1, int r1, int c1)
|
||||
{
|
||||
int rows = m1.rows();
|
||||
int cols = m1.cols();
|
||||
|
||||
Matrix<Scalar, Dynamic, Dynamic> mi = m1.minor(0,0).eval();
|
||||
VERIFY_IS_APPROX(mi, m1.block(1,1,rows-1,cols-1));
|
||||
mi = m1.minor(r1,c1);
|
||||
VERIFY_IS_APPROX(mi.transpose(), m1.transpose().minor(c1,r1));
|
||||
//check operator(), both constant and non-constant, on minor()
|
||||
m1.minor(r1,c1)(0,0) = m1.minor(0,0)(0,0);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Scalar> struct CheckMinor<Scalar,1,1>
|
||||
{
|
||||
typedef Matrix<Scalar, 1, 1> MatrixType;
|
||||
CheckMinor(MatrixType&, int, int) {}
|
||||
};
|
||||
|
||||
template<typename MatrixType> void block(const MatrixType& m)
|
||||
{
|
||||
/* this test covers the following files:
|
||||
Row.h Column.h Block.h Minor.h
|
||||
*/
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
|
||||
@ -101,9 +71,6 @@ template<typename MatrixType> void block(const MatrixType& m)
|
||||
m1.block(r1,c1,r2-r1+1,c2-c1+1) = s1 * m2.block(0, 0, r2-r1+1,c2-c1+1);
|
||||
m1.block(r1,c1,r2-r1+1,c2-c1+1)(r2-r1,c2-c1) = m2.block(0, 0, r2-r1+1,c2-c1+1)(0,0);
|
||||
|
||||
//check minor()
|
||||
CheckMinor<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> checkminor(m1,r1,c1);
|
||||
|
||||
const int BlockRows = EIGEN_ENUM_MIN(MatrixType::RowsAtCompileTime,2);
|
||||
const int BlockCols = EIGEN_ENUM_MIN(MatrixType::ColsAtCompileTime,5);
|
||||
if (rows>=5 && cols>=8)
|
||||
|
||||
@ -36,6 +36,7 @@ template<typename MatrixType> void inverse_permutation_4x4()
|
||||
MatrixType m = PermutationMatrix<4>(indices);
|
||||
MatrixType inv = m.inverse();
|
||||
double error = double( (m*inv-MatrixType::Identity()).norm() / NumTraits<Scalar>::epsilon() );
|
||||
EIGEN_DEBUG_VAR(error)
|
||||
VERIFY(error == 0.0);
|
||||
std::next_permutation(indices.data(),indices.data()+4);
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user