aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Core/Minor.h5
-rw-r--r--Eigen/src/LU/Determinant.h4
-rw-r--r--Eigen/src/LU/Inverse.h101
-rw-r--r--test/inverse.cpp11
-rw-r--r--test/main.h15
-rw-r--r--test/prec_inverse_4x4.cpp28
6 files changed, 43 insertions, 121 deletions
diff --git a/Eigen/src/Core/Minor.h b/Eigen/src/Core/Minor.h
index ab058b187..eb44f3760 100644
--- a/Eigen/src/Core/Minor.h
+++ b/Eigen/src/Core/Minor.h
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// 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
@@ -54,7 +54,8 @@ struct ei_traits<Minor<MatrixType> >
MaxColsAtCompileTime = (MatrixType::MaxColsAtCompileTime != Dynamic) ?
int(MatrixType::MaxColsAtCompileTime) - 1 : Dynamic,
Flags = _MatrixTypeNested::Flags & HereditaryBits,
- CoeffReadCost = _MatrixTypeNested::CoeffReadCost
+ CoeffReadCost = _MatrixTypeNested::CoeffReadCost // minor is used typically on tiny matrices,
+ // where loops are unrolled and the 'if' evaluates at compile time
};
};
diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
index 8870d9f20..27ad6abe9 100644
--- a/Eigen/src/LU/Determinant.h
+++ b/Eigen/src/LU/Determinant.h
@@ -118,7 +118,9 @@ template<typename Derived>
inline typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
{
assert(rows() == cols());
- return ei_determinant_impl<Derived>::run(derived());
+ typedef typename ei_nested<Derived,RowsAtCompileTime>::type Nested;
+ Nested nested(derived());
+ return ei_determinant_impl<typename ei_cleantype<Nested>::type>::run(nested);
}
#endif // EIGEN_DETERMINANT_H
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 9d5e86845..8afbfda96 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -183,92 +183,27 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
****************************/
template<typename MatrixType, typename ResultType>
-void ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType& result)
-{
- /* Let's split M into four 2x2 blocks:
- * (P Q)
- * (R S)
- * If P is invertible, with inverse denoted by P_inverse, and if
- * (S - R*P_inverse*Q) is also invertible, then the inverse of M is
- * (P' Q')
- * (R' S')
- * where
- * S' = (S - R*P_inverse*Q)^(-1)
- * P' = P1 + (P1*Q) * S' *(R*P_inverse)
- * Q' = -(P_inverse*Q) * S'
- * R' = -S' * (R*P_inverse)
- */
- typedef Block<ResultType,2,2> XprBlock22;
- typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22;
- Block22 P_inverse;
- ei_compute_inverse<XprBlock22, Block22>::run(matrix.template block<2,2>(0,0), P_inverse);
- const Block22 Q = matrix.template block<2,2>(0,2);
- const Block22 P_inverse_times_Q = P_inverse * Q;
- const XprBlock22 R = matrix.template block<2,2>(2,0);
- const Block22 R_times_P_inverse = R * P_inverse;
- const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q;
- const XprBlock22 S = matrix.template block<2,2>(2,2);
- const Block22 X = S - R_times_P_inverse_times_Q;
- Block22 Y;
- ei_compute_inverse<Block22, Block22>::run(X, Y);
- result.template block<2,2>(2,2) = Y;
- result.template block<2,2>(2,0) = - Y * R_times_P_inverse;
- const Block22 Z = P_inverse_times_Q * Y;
- result.template block<2,2>(0,2) = - Z;
- result.template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
-}
-
-template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 4>
{
- static inline void run(const MatrixType& _matrix, ResultType& result)
+ static inline void run(const MatrixType& matrix, ResultType& result)
{
- typedef typename ResultType::Scalar Scalar;
- typedef typename MatrixType::RealScalar RealScalar;
-
- // we will do row permutations on the matrix. This copy should have negligible cost.
- // if not, consider working in-place on the matrix (const-cast it, but then undo the permutations
- // to nevertheless honor constness)
- typename MatrixType::PlainMatrixType matrix(_matrix);
-
- // let's extract from the 2 first colums a 2x2 block whose determinant is as big as possible.
- int good_row0, good_row1, good_i;
- Matrix<RealScalar,6,1> absdet;
-
- // any 2x2 block with determinant above this threshold will be considered good enough.
- // The magic value 1e-1 here comes from experimentation. The bigger it is, the higher the precision,
- // the slower the computation. This value 1e-1 gives precision almost as good as the brutal cofactors
- // algorithm, both in average and in worst-case precision.
- RealScalar d = (matrix.col(0).squaredNorm()+matrix.col(1).squaredNorm()) * RealScalar(1e-1);
- #define ei_inv_size4_helper_macro(i,row0,row1) \
- absdet[i] = ei_abs(matrix.coeff(row0,0)*matrix.coeff(row1,1) \
- - matrix.coeff(row0,1)*matrix.coeff(row1,0)); \
- if(absdet[i] > d) { good_row0=row0; good_row1=row1; goto good; }
- ei_inv_size4_helper_macro(0,0,1)
- ei_inv_size4_helper_macro(1,0,2)
- ei_inv_size4_helper_macro(2,0,3)
- ei_inv_size4_helper_macro(3,1,2)
- ei_inv_size4_helper_macro(4,1,3)
- ei_inv_size4_helper_macro(5,2,3)
-
- // no 2x2 block has determinant bigger than the threshold. So just take the one that
- // has the biggest determinant
- absdet.maxCoeff(&good_i);
- good_row0 = good_i <= 2 ? 0 : good_i <= 4 ? 1 : 2;
- good_row1 = good_i <= 2 ? good_i+1 : good_i <= 4 ? good_i-1 : 3;
-
- // now good_row0 and good_row1 are correctly set
- good:
-
- // do row permutations to move this 2x2 block to the top
- matrix.row(0).swap(matrix.row(good_row0));
- matrix.row(1).swap(matrix.row(good_row1));
- // now applying our helper function is numerically stable
- ei_compute_inverse_size4_helper(matrix, result);
- // Since we did row permutations on the original matrix, we need to do column permutations
- // in the reverse order on the inverse
- result.col(1).swap(result.col(good_row1));
- result.col(0).swap(result.col(good_row0));
+ 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 /= (matrix.col(0).cwise()*result.row(0).transpose()).sum();
}
};
diff --git a/test/inverse.cpp b/test/inverse.cpp
index 59b791507..b80e139e0 100644
--- a/test/inverse.cpp
+++ b/test/inverse.cpp
@@ -38,18 +38,11 @@ template<typename MatrixType> void inverse(const MatrixType& m)
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
- MatrixType m1 = MatrixType::Random(rows, cols),
+ MatrixType m1(rows, cols),
m2(rows, cols),
mzero = MatrixType::Zero(rows, cols),
identity = MatrixType::Identity(rows, rows);
-
- if (ei_is_same_type<RealScalar,float>::ret)
- {
- // let's build a more stable to inverse matrix
- MatrixType a = MatrixType::Random(rows,cols);
- m1 += m1 * m1.adjoint() + a * a.adjoint();
- }
-
+ createRandomMatrixOfRank(rows,rows,rows,m1);
m2 = m1.inverse();
VERIFY_IS_APPROX(m1, m2.inverse() );
diff --git a/test/main.h b/test/main.h
index 8bc71b3ae..06c17a9ae 100644
--- a/test/main.h
+++ b/test/main.h
@@ -353,13 +353,26 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType&
typedef Matrix<Scalar, Rows, Rows> MatrixAType;
typedef Matrix<Scalar, Cols, Cols> MatrixBType;
+ if(desired_rank == 0)
+ {
+ m.setZero(rows,cols);
+ return;
+ }
+
+ if(desired_rank == 1)
+ {
+ m = VectorType::Random(rows) * VectorType::Random(cols).transpose();
+ return;
+ }
+
MatrixAType a = MatrixAType::Random(rows,rows);
MatrixType d = MatrixType::Identity(rows,cols);
MatrixBType b = MatrixBType::Random(cols,cols);
// set the diagonal such that only desired_rank non-zero entries reamain
const int diag_size = std::min(d.rows(),d.cols());
- d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank);
+ if(diag_size != desired_rank)
+ d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank);
HouseholderQR<MatrixAType> qra(a);
HouseholderQR<MatrixBType> qrb(b);
diff --git a/test/prec_inverse_4x4.cpp b/test/prec_inverse_4x4.cpp
index 613535346..83b1a8a71 100644
--- a/test/prec_inverse_4x4.cpp
+++ b/test/prec_inverse_4x4.cpp
@@ -26,28 +26,6 @@
#include <Eigen/LU>
#include <algorithm>
-Matrix4f inverse(const Matrix4f& m)
-{
- Matrix4f r;
- r(0,0) = m.minor(0,0).determinant();
- r(1,0) = -m.minor(0,1).determinant();
- r(2,0) = m.minor(0,2).determinant();
- r(3,0) = -m.minor(0,3).determinant();
- r(0,2) = m.minor(2,0).determinant();
- r(1,2) = -m.minor(2,1).determinant();
- r(2,2) = m.minor(2,2).determinant();
- r(3,2) = -m.minor(2,3).determinant();
- r(0,1) = -m.minor(1,0).determinant();
- r(1,1) = m.minor(1,1).determinant();
- r(2,1) = -m.minor(1,2).determinant();
- r(3,1) = m.minor(1,3).determinant();
- r(0,3) = -m.minor(3,0).determinant();
- r(1,3) = m.minor(3,1).determinant();
- r(2,3) = -m.minor(3,2).determinant();
- r(3,3) = m.minor(3,3).determinant();
- return r / (m(0,0)*r(0,0) + m(1,0)*r(0,1) + m(2,0)*r(0,2) + m(3,0)*r(0,3));
-}
-
template<typename MatrixType> void inverse_permutation_4x4()
{
typedef typename MatrixType::Scalar Scalar;
@@ -79,7 +57,7 @@ template<typename MatrixType> void inverse_general_4x4(int repeat)
do {
m = MatrixType::Random();
absdet = ei_abs(m.determinant());
- } while(absdet < 10 * epsilon<Scalar>());
+ } while(absdet < epsilon<Scalar>());
MatrixType inv = m.inverse();
double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() );
error_sum += error;
@@ -89,8 +67,8 @@ template<typename MatrixType> void inverse_general_4x4(int repeat)
double error_avg = error_sum / repeat;
EIGEN_DEBUG_VAR(error_avg);
EIGEN_DEBUG_VAR(error_max);
- VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.4 : 1.4) );
- VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 160.0 : 75.) );
+ VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.0));
+ VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 16.0));
}
void test_prec_inverse_4x4()
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-01 05:55:37 +0000