diff options
author | 2009-12-09 12:43:25 -0500 | |
---|---|---|
committer | 2009-12-09 12:43:25 -0500 | |
commit | d2e44f263631981d9e547caafe36b1de5ba785f9 (patch) | |
tree | 574d7aff0554739bea2cede8dc706cd5cd8c7c4a /Eigen/src/LU/Inverse.h | |
parent | f0315295e9ae2fd8afdc05d3e5b790b4660ffc58 (diff) |
* 4x4 inverse: revert to cofactors method
* inverse tests: use createRandomMatrixOfRank, use more strict precision
* tests: createRandomMatrixOfRank: support 1x1 matrices
* determinant: nest the xpr
* Minor: add comment
Diffstat (limited to 'Eigen/src/LU/Inverse.h')
-rw-r--r-- | Eigen/src/LU/Inverse.h | 101 |
1 files changed, 18 insertions, 83 deletions
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(); } }; |