From 5eabf2b75da885d8e0134e0c02f9b5f5f7dfa8c0 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 29 Jun 2009 00:00:29 +0200 Subject: [PATCH 01/25] double precision() : change to 1e-12 instead of 1e-11 (as discussed several times on the list) --- Eigen/src/Core/MathFunctions.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 8cb1f7b40..3527042f9 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -151,7 +151,7 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision inline double precision() { return 1e-11; } +template<> inline double precision() { return 1e-12; } template<> inline double machine_epsilon() { return 2.220e-16; } inline double ei_real(double x) { return x; } From bf91003d37f18d5bcf66bafb7290560f11cdf6ba Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 29 Jun 2009 00:08:34 +0200 Subject: [PATCH 02/25] FreeBSD: determine precisely when malloc is 16-byte aligned --- Eigen/src/Core/util/Memory.h | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index cc3aa4fac..f8581eebc 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -27,7 +27,17 @@ #ifndef EIGEN_MEMORY_H #define EIGEN_MEMORY_H -#if defined(__APPLE__) || defined(_WIN64) || defined (__FreeBSD__) +// FreeBSD 6 seems to have 16-byte aligned malloc +// See http://svn.freebsd.org/viewvc/base/stable/6/lib/libc/stdlib/malloc.c?view=markup +// FreeBSD 7 seems to have 16-byte aligned malloc except on ARM and MIPS architectures +// See http://svn.freebsd.org/viewvc/base/stable/7/lib/libc/stdlib/malloc.c?view=markup +#if defined(__FreeBSD__) && !defined(__arm__) && !defined(__mips__) +#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 1 +#else +#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 0 +#endif + +#if defined(__APPLE__) || defined(_WIN64) || EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED #define EIGEN_MALLOC_ALREADY_ALIGNED 1 #else #define EIGEN_MALLOC_ALREADY_ALIGNED 0 From 632cb7a4a1055fa3abbd341466d240cc50069361 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 29 Jun 2009 19:26:20 +0200 Subject: [PATCH 03/25] patch by Myguel from the forum: fix documentation --- Eigen/src/Geometry/Quaternion.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 66580cf1b..d7480bf9c 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -303,7 +303,9 @@ inline Quaternion& Quaternion::operator=(const MatrixBase inline typename Quaternion::Matrix3 Quaternion::toRotationMatrix(void) const From 126a031a391390f42d95f2adb9a13db6345ab007 Mon Sep 17 00:00:00 2001 From: Manuel Yguel Date: Mon, 29 Jun 2009 20:47:37 +0200 Subject: [PATCH 04/25] computeInverseWithCheck method added to matrix base (specialization for 1D to 4D) --- Eigen/src/Core/MatrixBase.h | 1 + Eigen/src/LU/Inverse.h | 201 +++++++++++++++++++++++++++++++++--- 2 files changed, 185 insertions(+), 17 deletions(-) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index addf59570..df411da31 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -643,6 +643,7 @@ template class MatrixBase const PartialLU partialLu() const; const PlainMatrixType inverse() const; void computeInverse(PlainMatrixType *result) const; + bool computeInverseWithCheck( PlainMatrixType *result ) const; Scalar determinant() const; /////////// Cholesky module /////////// diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index b29efb60c..02e23b407 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -29,15 +29,23 @@ *** Part 1 : optimized implementations for fixed-size 2,3,4 cases *** ********************************************************************/ +template +inline void ei_compute_inverse_in_size2_aux( + const XprType& matrix, const typename MatrixType::Scalar invdet, + MatrixType* result) +{ + result->coeffRef(0,0) = matrix.coeff(1,1) * invdet; + result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet; + result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet; + result->coeffRef(1,1) = matrix.coeff(0,0) * invdet; +} + template void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result) { typedef typename MatrixType::Scalar Scalar; const Scalar invdet = Scalar(1) / matrix.determinant(); - result->coeffRef(0,0) = matrix.coeff(1,1) * invdet; - result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet; - result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet; - result->coeffRef(1,1) = matrix.coeff(0,0) * invdet; + ei_compute_inverse_in_size2_aux( matrix, invdet, result ); } template @@ -47,13 +55,31 @@ bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixTy const Scalar det = matrix.determinant(); if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false; const Scalar invdet = Scalar(1) / det; - result->coeffRef(0,0) = matrix.coeff(1,1) * invdet; - result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet; - result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet; - result->coeffRef(1,1) = matrix.coeff(0,0) * invdet; + ei_compute_inverse_in_size2_aux( matrix, invdet, result ); return true; } +template +inline void ei_compute_inverse_in_size3_aux( + const XprType& matrix, + const typename MatrixType::Scalar invdet, + const typename MatrixType::Scalar det_minor00, + const typename MatrixType::Scalar det_minor10, + const typename MatrixType::Scalar det_minor20, + MatrixType* result) +{ + result->coeffRef(0, 0) = det_minor00 * invdet; + result->coeffRef(0, 1) = -det_minor10 * invdet; + result->coeffRef(0, 2) = det_minor20 * 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; +} + + template void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result) { @@ -64,15 +90,23 @@ void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* resu const Scalar invdet = Scalar(1) / ( det_minor00 * matrix.coeff(0,0) - det_minor10 * matrix.coeff(1,0) + det_minor20 * matrix.coeff(2,0) ); - result->coeffRef(0, 0) = det_minor00 * invdet; - result->coeffRef(0, 1) = -det_minor10 * invdet; - result->coeffRef(0, 2) = det_minor20 * 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; + ei_compute_inverse_in_size3_aux( matrix, invdet, det_minor00, det_minor10, det_minor20, result ); +} + +template +bool ei_compute_inverse_in_size3_case_with_check(const XprType& matrix, MatrixType* result) +{ + typedef typename MatrixType::Scalar Scalar; + 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(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false; + const Scalar invdet = Scalar(1) / det; + ei_compute_inverse_in_size3_aux( matrix, invdet, det_minor00, det_minor10, det_minor20, result ); + return true; } template @@ -161,6 +195,59 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu } } + +template +bool ei_compute_inverse_in_size4_case_with_check(const XprType& matrix, MatrixType* result) +{ + if(ei_compute_inverse_in_size4_case_helper(matrix, result)) + { + // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful. + return true; + } + else + { + // rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be). + // since this is a rare case, we don't need to optimize it. We just want to handle it with little + // additional code. + MatrixType m(matrix); + m.row(0).swap(m.row(2)); + m.row(1).swap(m.row(3)); + if(ei_compute_inverse_in_size4_case_helper(m, result)) + { + // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some + // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting + // the corresponding columns. + result->col(0).swap(result->col(2)); + result->col(1).swap(result->col(3)); + return true; + } + else + { + // last possible case. Since matrix is assumed to be invertible, this last case has to work. + // first, undo the swaps previously made + m.row(0).swap(m.row(2)); + m.row(1).swap(m.row(3)); + // swap row 0 with the the row among 0 and 1 that has the biggest 2 first coeffs + int swap0with = ei_abs(m.coeff(0,0))+ei_abs(m.coeff(0,1))>ei_abs(m.coeff(1,0))+ei_abs(m.coeff(1,1)) ? 0 : 1; + m.row(0).swap(m.row(swap0with)); + // swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs + int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3; + m.row(1).swap(m.row(swap1with)); + if( ei_compute_inverse_in_size4_case_helper(m, result) ) + { + result->col(1).swap(result->col(swap1with)); + result->col(0).swap(result->col(swap0with)); + return true; + } + else{ + return false; } + } + } + +} + + + /*********************************************** *** Part 2 : selector and MatrixBase methods *** ***********************************************/ @@ -254,4 +341,84 @@ inline const typename MatrixBase::PlainMatrixType MatrixBase:: return result; } + +/***************************************** + * Compute inverse with invertibility check +*********************************************/ + +template +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + typedef typename MatrixType::Scalar Scalar; + LU lu( matrix ); + if( !lu.isInvertible() ) return false; + lu.computeInverse(result); + return true; + } +}; + +template +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + if( 0 == result->coeffRef(0,0) ) return false; + + typedef typename MatrixType::Scalar Scalar; + result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0); + return true; + } +}; + +template +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_in_size2_case_with_check(matrix, result); + } +}; + +template +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_in_size3_case_with_check(matrix, result); + } +}; + +template +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_in_size4_case_with_check(matrix, result); + } +}; + +/** \lu_module + * + * If the matrix is invertible, computes the matrix inverse of this matrix + * and returns true otherwise return false. + * + * \note This matrix must be invertible, otherwise the result is undefined. + * + * \param result Pointer to the matrix in which to store the result. Undefined + * if the matrix is not invertible. + * \return true if the matrix is invertible false otherwise. + * + * \sa inverse() + */ +template +inline bool MatrixBase::computeInverseWithCheck(PlainMatrixType *result) const +{ + ei_assert(rows() == cols()); + EIGEN_STATIC_ASSERT(NumTraits::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) + return ei_compute_inverse_with_check::run(eval(), result); +} + + #endif // EIGEN_INVERSE_H From 7b750182f244683d37de585b6181cf59bf24ff2b Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 29 Jun 2009 22:07:37 +0200 Subject: [PATCH 05/25] * polish computeInverseWithCheck to share more code, fix documentation, fix coding style * add snippet for computeInverseWithCheck documentation * expand unit-tests to cover computeInverseWithCheck --- Eigen/src/Core/MatrixBase.h | 2 +- Eigen/src/LU/Inverse.h | 200 +++++++----------- .../MatrixBase_computeInverseWithCheck.cpp | 9 + test/inverse.cpp | 4 + test/lu.cpp | 5 + 5 files changed, 90 insertions(+), 130 deletions(-) create mode 100644 doc/snippets/MatrixBase_computeInverseWithCheck.cpp diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index df411da31..f925e0b0f 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -643,7 +643,7 @@ template class MatrixBase const PartialLU partialLu() const; const PlainMatrixType inverse() const; void computeInverse(PlainMatrixType *result) const; - bool computeInverseWithCheck( PlainMatrixType *result ) const; + bool computeInverseWithCheck( PlainMatrixType *result ) const; Scalar determinant() const; /////////// Cholesky module /////////// diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 02e23b407..5af14813d 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -30,43 +30,43 @@ ********************************************************************/ template -inline void ei_compute_inverse_in_size2_aux( - const XprType& matrix, const typename MatrixType::Scalar invdet, - MatrixType* result) +inline void ei_compute_inverse_size2_helper( + const XprType& matrix, const typename MatrixType::Scalar& invdet, + MatrixType* result) { - result->coeffRef(0,0) = matrix.coeff(1,1) * invdet; - result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet; - result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet; - result->coeffRef(1,1) = matrix.coeff(0,0) * invdet; + result->coeffRef(0,0) = matrix.coeff(1,1) * invdet; + result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet; + result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet; + result->coeffRef(1,1) = matrix.coeff(0,0) * invdet; } template -void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result) +inline void ei_compute_inverse_size2(const MatrixType& matrix, MatrixType* result) { typedef typename MatrixType::Scalar Scalar; const Scalar invdet = Scalar(1) / matrix.determinant(); - ei_compute_inverse_in_size2_aux( matrix, invdet, result ); + ei_compute_inverse_size2_helper( matrix, invdet, result ); } template -bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixType* result) +bool ei_compute_inverse_size2_with_check(const XprType& matrix, MatrixType* result) { typedef typename MatrixType::Scalar Scalar; const Scalar det = matrix.determinant(); if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false; const Scalar invdet = Scalar(1) / det; - ei_compute_inverse_in_size2_aux( matrix, invdet, result ); + ei_compute_inverse_size2_helper( matrix, invdet, result ); return true; } template -inline void ei_compute_inverse_in_size3_aux( - const XprType& matrix, - const typename MatrixType::Scalar invdet, - const typename MatrixType::Scalar det_minor00, - const typename MatrixType::Scalar det_minor10, - const typename MatrixType::Scalar det_minor20, - MatrixType* result) +void ei_compute_inverse_size3_helper( + const XprType& matrix, + const typename MatrixType::Scalar& invdet, + const typename MatrixType::Scalar& det_minor00, + const typename MatrixType::Scalar& det_minor10, + const typename MatrixType::Scalar& det_minor20, + MatrixType* result) { result->coeffRef(0, 0) = det_minor00 * invdet; result->coeffRef(0, 1) = -det_minor10 * invdet; @@ -79,38 +79,24 @@ inline void ei_compute_inverse_in_size3_aux( result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet; } - -template -void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result) -{ - typedef typename MatrixType::Scalar Scalar; - 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 invdet = Scalar(1) / ( det_minor00 * matrix.coeff(0,0) - - det_minor10 * matrix.coeff(1,0) - + det_minor20 * matrix.coeff(2,0) ); - ei_compute_inverse_in_size3_aux( matrix, invdet, det_minor00, det_minor10, det_minor20, result ); -} - -template -bool ei_compute_inverse_in_size3_case_with_check(const XprType& matrix, MatrixType* result) +template +bool ei_compute_inverse_size3(const XprType& matrix, MatrixType* result) { typedef typename MatrixType::Scalar Scalar; 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(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false; + - det_minor10 * matrix.coeff(1,0) + + det_minor20 * matrix.coeff(2,0) ); + if(Check) if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false; const Scalar invdet = Scalar(1) / det; - ei_compute_inverse_in_size3_aux( matrix, invdet, det_minor00, det_minor10, det_minor20, result ); + ei_compute_inverse_size3_helper( matrix, invdet, det_minor00, det_minor10, det_minor20, result ); return true; } template -bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result) +bool ei_compute_inverse_size4_helper(const MatrixType& matrix, MatrixType* result) { /* Let's split M into four 2x2 blocks: * (P Q) @@ -128,7 +114,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp typedef Block XprBlock22; typedef typename MatrixBase::PlainMatrixType Block22; Block22 P_inverse; - if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse)) + if(ei_compute_inverse_size2_with_check(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; @@ -138,7 +124,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp 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_in_size2_case(X, &Y); + ei_compute_inverse_size2(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; @@ -152,54 +138,10 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp } } -template -void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result) -{ - if(ei_compute_inverse_in_size4_case_helper(matrix, result)) - { - // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful. - return; - } - else - { - // rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be). - // since this is a rare case, we don't need to optimize it. We just want to handle it with little - // additional code. - MatrixType m(matrix); - m.row(0).swap(m.row(2)); - m.row(1).swap(m.row(3)); - if(ei_compute_inverse_in_size4_case_helper(m, result)) - { - // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some - // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting - // the corresponding columns. - result->col(0).swap(result->col(2)); - result->col(1).swap(result->col(3)); - } - else - { - // last possible case. Since matrix is assumed to be invertible, this last case has to work. - // first, undo the swaps previously made - m.row(0).swap(m.row(2)); - m.row(1).swap(m.row(3)); - // swap row 0 with the the row among 0 and 1 that has the biggest 2 first coeffs - int swap0with = ei_abs(m.coeff(0,0))+ei_abs(m.coeff(0,1))>ei_abs(m.coeff(1,0))+ei_abs(m.coeff(1,1)) ? 0 : 1; - m.row(0).swap(m.row(swap0with)); - // swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs - int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3; - m.row(1).swap(m.row(swap1with)); - ei_compute_inverse_in_size4_case_helper(m, result); - result->col(1).swap(result->col(swap1with)); - result->col(0).swap(result->col(swap0with)); - } - } -} - - template -bool ei_compute_inverse_in_size4_case_with_check(const XprType& matrix, MatrixType* result) +bool ei_compute_inverse_size4_with_check(const XprType& matrix, MatrixType* result) { - if(ei_compute_inverse_in_size4_case_helper(matrix, result)) + if(ei_compute_inverse_size4_helper(matrix, result)) { // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful. return true; @@ -212,18 +154,17 @@ bool ei_compute_inverse_in_size4_case_with_check(const XprType& matrix, MatrixTy MatrixType m(matrix); m.row(0).swap(m.row(2)); m.row(1).swap(m.row(3)); - if(ei_compute_inverse_in_size4_case_helper(m, result)) + if(ei_compute_inverse_size4_helper(m, result)) { // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting // the corresponding columns. result->col(0).swap(result->col(2)); result->col(1).swap(result->col(3)); - return true; + return true; } else { - // last possible case. Since matrix is assumed to be invertible, this last case has to work. // first, undo the swaps previously made m.row(0).swap(m.row(2)); m.row(1).swap(m.row(3)); @@ -233,17 +174,19 @@ bool ei_compute_inverse_in_size4_case_with_check(const XprType& matrix, MatrixTy // swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3; m.row(1).swap(m.row(swap1with)); - if( ei_compute_inverse_in_size4_case_helper(m, result) ) - { - result->col(1).swap(result->col(swap1with)); - result->col(0).swap(result->col(swap0with)); - return true; - } - else{ - return false; } + if( ei_compute_inverse_size4_helper(m, result) ) + { + result->col(1).swap(result->col(swap1with)); + result->col(0).swap(result->col(swap0with)); + return true; + } + else + { + // non-invertible matrix + return false; + } } } - } @@ -276,7 +219,7 @@ struct ei_compute_inverse { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_size2_case(matrix, result); + ei_compute_inverse_size2(matrix, result); } }; @@ -285,7 +228,7 @@ struct ei_compute_inverse { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_size3_case(matrix, result); + ei_compute_inverse_size3(matrix, result); } }; @@ -294,7 +237,7 @@ struct ei_compute_inverse { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_size4_case(matrix, result); + ei_compute_inverse_size4_with_check(matrix, result); } }; @@ -302,14 +245,15 @@ struct ei_compute_inverse * * Computes the matrix inverse of this matrix. * - * \note This matrix must be invertible, otherwise the result is undefined. + * \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use + * computeInverseWithCheck(). * * \param result Pointer to the matrix in which to store the result. * * Example: \include MatrixBase_computeInverse.cpp * Output: \verbinclude MatrixBase_computeInverse.out * - * \sa inverse() + * \sa inverse(), computeInverseWithCheck() */ template inline void MatrixBase::computeInverse(PlainMatrixType *result) const @@ -323,7 +267,8 @@ inline void MatrixBase::computeInverse(PlainMatrixType *result) const * * \returns the matrix inverse of this matrix. * - * \note This matrix must be invertible, otherwise the result is undefined. + * \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use + * computeInverseWithCheck(). * * \note This method returns a matrix by value, which can be inefficient. To avoid that overhead, * use computeInverse() instead. @@ -331,7 +276,7 @@ inline void MatrixBase::computeInverse(PlainMatrixType *result) const * Example: \include MatrixBase_inverse.cpp * Output: \verbinclude MatrixBase_inverse.out * - * \sa computeInverse() + * \sa computeInverse(), computeInverseWithCheck() */ template inline const typename MatrixBase::PlainMatrixType MatrixBase::inverse() const @@ -342,20 +287,20 @@ inline const typename MatrixBase::PlainMatrixType MatrixBase:: } -/***************************************** - * Compute inverse with invertibility check -*********************************************/ +/******************************************** + * Compute inverse with invertibility check * + *******************************************/ template struct ei_compute_inverse_with_check { static inline bool run(const MatrixType& matrix, MatrixType* result) { - typedef typename MatrixType::Scalar Scalar; - LU lu( matrix ); - if( !lu.isInvertible() ) return false; - lu.computeInverse(result); - return true; + typedef typename MatrixType::Scalar Scalar; + LU lu( matrix ); + if( !lu.isInvertible() ) return false; + lu.computeInverse(result); + return true; } }; @@ -364,11 +309,11 @@ struct ei_compute_inverse_with_check { static inline bool run(const MatrixType& matrix, MatrixType* result) { - if( 0 == result->coeffRef(0,0) ) return false; - - typedef typename MatrixType::Scalar Scalar; - result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0); - return true; + if( 0 == result->coeffRef(0,0) ) return false; + + typedef typename MatrixType::Scalar Scalar; + result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0); + return true; } }; @@ -377,7 +322,7 @@ struct ei_compute_inverse_with_check { static inline bool run(const MatrixType& matrix, MatrixType* result) { - return ei_compute_inverse_in_size2_case_with_check(matrix, result); + return ei_compute_inverse_size2_with_check(matrix, result); } }; @@ -386,7 +331,7 @@ struct ei_compute_inverse_with_check { static inline bool run(const MatrixType& matrix, MatrixType* result) { - return ei_compute_inverse_in_size3_case_with_check(matrix, result); + return ei_compute_inverse_size3(matrix, result); } }; @@ -395,22 +340,19 @@ struct ei_compute_inverse_with_check { static inline bool run(const MatrixType& matrix, MatrixType* result) { - return ei_compute_inverse_in_size4_case_with_check(matrix, result); + return ei_compute_inverse_size4_with_check(matrix, result); } }; /** \lu_module * - * If the matrix is invertible, computes the matrix inverse of this matrix - * and returns true otherwise return false. + * Computation of matrix inverse, with invertibility check. * - * \note This matrix must be invertible, otherwise the result is undefined. + * \returns true if the matrix is invertible, false otherwise. * - * \param result Pointer to the matrix in which to store the result. Undefined - * if the matrix is not invertible. - * \return true if the matrix is invertible false otherwise. + * \param result Pointer to the matrix in which to store the result. * - * \sa inverse() + * \sa inverse(), computeInverse() */ template inline bool MatrixBase::computeInverseWithCheck(PlainMatrixType *result) const diff --git a/doc/snippets/MatrixBase_computeInverseWithCheck.cpp b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp new file mode 100644 index 000000000..19e24c90b --- /dev/null +++ b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp @@ -0,0 +1,9 @@ +Matrix3d m = Matrix3d::Random(); +cout << "Here is the matrix m:" << endl << m << endl; +Matrix3d inv; +if(m.computeInverseWithCheck(&inv)) { + cout << "It is invertible, and its inverse is:" << endl << inv << endl; +} +else { + cout << "It is not invertible." << endl; +} diff --git a/test/inverse.cpp b/test/inverse.cpp index 90b897197..cdac7cdec 100644 --- a/test/inverse.cpp +++ b/test/inverse.cpp @@ -65,6 +65,10 @@ template void inverse(const MatrixType& m) // since for the general case we implement separately row-major and col-major, test that VERIFY_IS_APPROX(m1.transpose().inverse(), m1.inverse().transpose()); + + bool invertible = m1.computeInverseWithCheck(&m2); + VERIFY(invertible); + VERIFY_IS_APPROX(identity, m1*m2); } void test_inverse() diff --git a/test/lu.cpp b/test/lu.cpp index 60e45ff2b..4ad92bb11 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -57,6 +57,11 @@ template void lu_non_invertible() VERIFY_IS_APPROX(m3, m1*m2); m3 = MatrixType::Random(rows,cols2); VERIFY(!lu.solve(m3, &m2)); + + typedef Matrix SquareMatrixType; + SquareMatrixType m4(rows, rows), m5(rows, rows); + createRandomMatrixOfRank(rows/2, rows, rows, m4); + VERIFY(!m4.computeInverseWithCheck(&m5)); } template void lu_invertible() From d457ec58108ecfdaa8a8a5cedb2600828019077e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 4 Jul 2009 09:28:11 +0200 Subject: [PATCH 06/25] fix #20: SVD::solve() now resize the result --- Eigen/src/SVD/SVD.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 9733fbd21..45a7fbfa5 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -554,6 +554,8 @@ bool SVD::solve(const MatrixBase &b, ResultType* resul const int rows = m_matU.rows(); ei_assert(b.rows() == rows); + result->resize(m_matV.rows(), b.cols()); + Scalar maxVal = m_sigma.cwise().abs().maxCoeff(); for (int j=0; j Date: Sat, 4 Jul 2009 14:55:25 +0200 Subject: [PATCH 07/25] another test in the non invertible case --- test/inverse.cpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/test/inverse.cpp b/test/inverse.cpp index cdac7cdec..5ac39e35a 100644 --- a/test/inverse.cpp +++ b/test/inverse.cpp @@ -65,10 +65,21 @@ template void inverse(const MatrixType& m) // since for the general case we implement separately row-major and col-major, test that VERIFY_IS_APPROX(m1.transpose().inverse(), m1.inverse().transpose()); - + + //computeInverseWithCheck tests + //First: an invertible matrix bool invertible = m1.computeInverseWithCheck(&m2); VERIFY(invertible); VERIFY_IS_APPROX(identity, m1*m2); + + //Second: a rank one matrix (not invertible, except for 1x1 matrices) + VectorType v3 = VectorType::Random(rows); + MatrixType m3 = v3*v3.transpose(), m4(rows,cols); + invertible = m3.computeInverseWithCheck( &m4 ); + if( 1 == rows ){ + VERIFY( invertible ); } + else{ + VERIFY( !invertible ); } } void test_inverse() @@ -81,7 +92,7 @@ void test_inverse() CALL_SUBTEST( inverse(MatrixXf(8,8)) ); CALL_SUBTEST( inverse(MatrixXcd(7,7)) ); } - + // test some tricky cases for 4x4 matrices VERIFY_IS_APPROX((Matrix4f() << 0,0,1,0, 1,0,0,0, 0,1,0,0, 0,0,0,1).finished().inverse(), (Matrix4f() << 0,1,0,0, 0,0,1,0, 1,0,0,0, 0,0,0,1).finished()); From 60467e54a5650b5836b5f35bab9730527d250b70 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 5 Jul 2009 01:52:42 +0200 Subject: [PATCH 08/25] some docs improvements --- doc/D09_StructHavingEigenMembers.dox | 18 +++--------------- doc/D11_UnalignedArrayAssert.dox | 12 ++++++++++++ 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/doc/D09_StructHavingEigenMembers.dox b/doc/D09_StructHavingEigenMembers.dox index 4ff2b23aa..07b6570dc 100644 --- a/doc/D09_StructHavingEigenMembers.dox +++ b/doc/D09_StructHavingEigenMembers.dox @@ -88,7 +88,7 @@ The solution is to let class Foo have an aligned "operator new", as we showed in \section movetotop Should I then put all the members of Eigen types at the beginning of my class? -No, that's not needed. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So when you have code like +That's not required. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So code like this works fine: \code class Foo @@ -100,25 +100,13 @@ public: }; \endcode -it will work just fine. You do \b not need to rewrite it as - -\code -class Foo -{ - Eigen::Vector2d v; - double x; -public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW -}; -\endcode - \section dynamicsize What about dynamic-size matrices and vectors? -Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with fixed-size matrices and vectors. +Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with \ref FixedSizeVectorizable "fixed-size vectorizable matrices and vectors". \section bugineigen So is this a bug in Eigen? -No, it's not our bug. It's more like an inherent problem of the C++ language -- though it must be said that any other existing language probably has the same problem. The problem is that there is no way that you can specify an aligned "operator new" that would propagate to classes having you as member data. +No, it's not our bug. It's more like an inherent problem of the C++98 language specification, and seems to be taken care of in the upcoming language revision: see this document. \section conditional What if I want to do this conditionnally (depending on template parameters) ? diff --git a/doc/D11_UnalignedArrayAssert.dox b/doc/D11_UnalignedArrayAssert.dox index da1b49a08..e95c767b4 100644 --- a/doc/D11_UnalignedArrayAssert.dox +++ b/doc/D11_UnalignedArrayAssert.dox @@ -15,12 +15,24 @@ is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html There are 4 known causes for this issue. Please read on to understand them and learn how to fix them. \b Table \b of \b contents + - \ref where - \ref c1 - \ref c2 - \ref c3 - \ref c4 - \ref explanation +\section where Where in my own code is the cause of the problem? + +First of all, you need to find out where in your own code this assertion was triggered from. At first glance, the error message doesn't look helpful, as it refers to a file inside Eigen! However, since your program crashed, if you can reproduce the crash, you can get a backtrace using any debugger. For example, if you're using GCC, you can use the GDB debugger as follows: +\code +$ gdb ./my_program # Start GDB on your program +> run # Start running your program +... # Now reproduce the crash! +> bt # Obtain the backtrace +\endcode +Now that you know precisely where in your own code the problem is happening, read on to understand what you need to change. + \section c1 Cause 1: Structures having Eigen objects as members If you have code like this, From eec334c604d3fbd0ab7fcb0d9380f87c943afc90 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 5 Jul 2009 10:48:57 +0200 Subject: [PATCH 09/25] fixes a segfault --- Eigen/src/Sparse/SparseMatrix.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 03bc60bb5..b751b5cb9 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -371,7 +371,7 @@ class SparseMatrix const int outerSize = IsRowMajor ? rows : cols; m_innerSize = IsRowMajor ? cols : rows; m_data.clear(); - if (m_outerSize != outerSize) + if (m_outerSize != outerSize || m_outerSize==0) { delete[] m_outerIndex; m_outerIndex = new int [outerSize+1]; From c6f610093bff10c5a6e4fac94bfed422ba54b39a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 5 Jul 2009 11:33:55 +0200 Subject: [PATCH 10/25] add a VectorBlock expr as a specialization of Block --- Eigen/Core | 1 + Eigen/src/Core/Block.h | 243 ----------------- Eigen/src/Core/MatrixBase.h | 24 +- Eigen/src/Core/VectorBlock.h | 311 ++++++++++++++++++++++ Eigen/src/Core/util/ForwardDeclarations.h | 1 + Eigen/src/Core/util/XprHelper.h | 2 - doc/examples/class_FixedVectorBlock.cpp | 26 ++ doc/examples/class_VectorBlock.cpp | 26 ++ 8 files changed, 377 insertions(+), 257 deletions(-) create mode 100644 Eigen/src/Core/VectorBlock.h create mode 100644 doc/examples/class_FixedVectorBlock.cpp create mode 100644 doc/examples/class_VectorBlock.cpp diff --git a/Eigen/Core b/Eigen/Core index ace3a204e..0d3650c99 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -164,6 +164,7 @@ namespace Eigen { #include "src/Core/MapBase.h" #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" diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h index de4268344..d7af971f9 100644 --- a/Eigen/src/Core/Block.h +++ b/Eigen/src/Core/Block.h @@ -314,249 +314,6 @@ inline const typename BlockReturnType::Type MatrixBase return typename BlockReturnType::Type(derived(), startRow, startCol, blockRows, blockCols); } -/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this. - * - * \only_for_vectors - * - * \addexample SegmentIntInt \label How to reference a sub-vector (dynamic size) - * - * \param start the first coefficient in the segment - * \param size the number of coefficients in the segment - * - * Example: \include MatrixBase_segment_int_int.cpp - * Output: \verbinclude MatrixBase_segment_int_int.out - * - * \note Even though the returned expression has dynamic size, in the case - * when it is applied to a fixed-size vector, it inherits a fixed maximal size, - * which means that evaluating it does not cause a dynamic memory allocation. - * - * \sa class Block, segment(int) - */ -template -inline typename BlockReturnType::SubVectorType MatrixBase - ::segment(int start, int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename BlockReturnType::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** This is the const version of segment(int,int).*/ -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::segment(int start, int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename BlockReturnType::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** \returns a dynamic-size expression of the first coefficients of *this. - * - * \only_for_vectors - * - * \param size the number of coefficients in the block - * - * \addexample BlockInt \label How to reference a sub-vector (fixed-size) - * - * Example: \include MatrixBase_start_int.cpp - * Output: \verbinclude MatrixBase_start_int.out - * - * \note Even though the returned expression has dynamic size, in the case - * when it is applied to a fixed-size vector, it inherits a fixed maximal size, - * which means that evaluating it does not cause a dynamic memory allocation. - * - * \sa class Block, block(int,int) - */ -template -inline typename BlockReturnType::SubVectorType -MatrixBase::start(int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), 0, 0, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** This is the const version of start(int).*/ -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::start(int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), 0, 0, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** \returns a dynamic-size expression of the last coefficients of *this. - * - * \only_for_vectors - * - * \param size the number of coefficients in the block - * - * \addexample BlockEnd \label How to reference the end of a vector (fixed-size) - * - * Example: \include MatrixBase_end_int.cpp - * Output: \verbinclude MatrixBase_end_int.out - * - * \note Even though the returned expression has dynamic size, in the case - * when it is applied to a fixed-size vector, it inherits a fixed maximal size, - * which means that evaluating it does not cause a dynamic memory allocation. - * - * \sa class Block, block(int,int) - */ -template -inline typename BlockReturnType::SubVectorType -MatrixBase::end(int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - size, - ColsAtCompileTime == 1 ? 0 : cols() - size, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** This is the const version of end(int).*/ -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::end(int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - size, - ColsAtCompileTime == 1 ? 0 : cols() - size, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this - * - * \only_for_vectors - * - * The template parameter \a Size is the number of coefficients in the block - * - * \param start the index of the first element of the sub-vector - * - * Example: \include MatrixBase_template_int_segment.cpp - * Output: \verbinclude MatrixBase_template_int_segment.out - * - * \sa class Block - */ -template -template -inline typename BlockReturnType::SubVectorType -MatrixBase::segment(int start) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start); -} - -/** This is the const version of segment(int).*/ -template -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::segment(int start) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start); -} - -/** \returns a fixed-size expression of the first coefficients of *this. - * - * \only_for_vectors - * - * The template parameter \a Size is the number of coefficients in the block - * - * \addexample BlockStart \label How to reference the start of a vector (fixed-size) - * - * Example: \include MatrixBase_template_int_start.cpp - * Output: \verbinclude MatrixBase_template_int_start.out - * - * \sa class Block - */ -template -template -inline typename BlockReturnType::SubVectorType -MatrixBase::start() -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block(derived(), 0, 0); -} - -/** This is the const version of start().*/ -template -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::start() const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block(derived(), 0, 0); -} - -/** \returns a fixed-size expression of the last coefficients of *this. - * - * \only_for_vectors - * - * The template parameter \a Size is the number of coefficients in the block - * - * Example: \include MatrixBase_template_int_end.cpp - * Output: \verbinclude MatrixBase_template_int_end.out - * - * \sa class Block - */ -template -template -inline typename BlockReturnType::SubVectorType -MatrixBase::end() -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - Size, - ColsAtCompileTime == 1 ? 0 : cols() - Size); -} - -/** This is the const version of end.*/ -template -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::end() const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - Size, - ColsAtCompileTime == 1 ? 0 : cols() - Size); -} - /** \returns a dynamic-size expression of a corner of *this. * * \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight, diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index f925e0b0f..65ab02d62 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -401,14 +401,14 @@ template class MatrixBase const typename BlockReturnType::Type block(int startRow, int startCol, int blockRows, int blockCols) const; - typename BlockReturnType::SubVectorType segment(int start, int size); - const typename BlockReturnType::SubVectorType segment(int start, int size) const; + VectorBlock segment(int start, int size); + const VectorBlock segment(int start, int size) const; - typename BlockReturnType::SubVectorType start(int size); - const typename BlockReturnType::SubVectorType start(int size) const; + VectorBlock start(int size); + const VectorBlock start(int size) const; - typename BlockReturnType::SubVectorType end(int size); - const typename BlockReturnType::SubVectorType end(int size) const; + VectorBlock end(int size); + const VectorBlock end(int size) const; typename BlockReturnType::Type corner(CornerType type, int cRows, int cCols); const typename BlockReturnType::Type corner(CornerType type, int cRows, int cCols) const; @@ -423,14 +423,14 @@ template class MatrixBase template const typename BlockReturnType::Type corner(CornerType type) const; - template typename BlockReturnType::SubVectorType start(void); - template const typename BlockReturnType::SubVectorType start() const; + template VectorBlock start(void); + template const VectorBlock start() const; - template typename BlockReturnType::SubVectorType end(); - template const typename BlockReturnType::SubVectorType end() const; + template VectorBlock end(); + template const VectorBlock end() const; - template typename BlockReturnType::SubVectorType segment(int start); - template const typename BlockReturnType::SubVectorType segment(int start) const; + template VectorBlock segment(int start); + template const VectorBlock segment(int start) const; Diagonal diagonal(); const Diagonal diagonal() const; diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h new file mode 100644 index 000000000..7ce5977f6 --- /dev/null +++ b/Eigen/src/Core/VectorBlock.h @@ -0,0 +1,311 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2006-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_VECTORBLOCK_H +#define EIGEN_VECTORBLOCK_H + +/** \class VectorBlock + * + * \brief Expression of a fixed-size or dynamic-size sub-vector + * + * \param VectorType the type of the object in which we are taking a sub-vector + * \param Size size of the sub-vector we are taking at compile time (optional) + * \param _PacketAccess allows to enforce aligned loads and stores if set to ForceAligned. + * The default is AsRequested. This parameter is internaly used by Eigen + * in expressions such as \code mat.segment() += other; \endcode and most of + * the time this is the only way it is used. + * + * This class represents an expression of either a fixed-size or dynamic-size sub-vector. + * It is the return type of MatrixBase::segment(int,int) and MatrixBase::segment(int) and + * most of the time this is the only way it is used. + * + * However, if you want to directly maniputate sub-vector expressions, + * for instance if you want to write a function returning such an expression, you + * will need to use this class. + * + * Here is an example illustrating the dynamic case: + * \include class_VectorBlock.cpp + * Output: \verbinclude class_VectorBlock.out + * + * \note Even though this expression has dynamic size, in the case where \a VectorType + * has fixed size, this expression inherits a fixed maximal size which means that evaluating + * it does not cause a dynamic memory allocation. + * + * Here is an example illustrating the fixed-size case: + * \include class_FixedVectorBlock.cpp + * Output: \verbinclude class_FixedVectorBlock.out + * + * \sa class Block, MatrixBase::segment(int,int,int,int), MatrixBase::segment(int,int) + */ +template +struct ei_traits > + : public ei_traits::RowsAtCompileTime==1 ? 1 : Size, + ei_traits::ColsAtCompileTime==1 ? 1 : Size, + _PacketAccess> > +{ +}; + +template class VectorBlock + : public Block::RowsAtCompileTime==1 ? 1 : Size, + ei_traits::ColsAtCompileTime==1 ? 1 : Size, + PacketAccess> +{ + typedef Block::RowsAtCompileTime==1 ? 1 : Size, + ei_traits::ColsAtCompileTime==1 ? 1 : Size, + PacketAccess> Base; + enum { + IsColVector = ei_traits::ColsAtCompileTime==1 + }; + public: + + using Base::operator=; + using Base::operator+=; + using Base::operator-=; + using Base::operator*=; + using Base::operator/=; + + /** Dynamic-size constructor + */ + inline VectorBlock(const VectorType& vector, int start, int size) + : Base(vector, + IsColVector ? start : 0, IsColVector ? 0 : start, + IsColVector ? size : 1, IsColVector ? 1 : size) + { + + EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock); + } + + /** Fixed-size constructor + */ + inline VectorBlock(const VectorType& vector, int start) + : Base(vector, IsColVector ? start : 0, IsColVector ? 0 : start) + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock); + } +}; + + +/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this. + * + * \only_for_vectors + * + * \addexample VectorBlockIntInt \label How to reference a sub-vector (dynamic size) + * + * \param start the first coefficient in the segment + * \param size the number of coefficients in the segment + * + * Example: \include MatrixBase_segment_int_int.cpp + * Output: \verbinclude MatrixBase_segment_int_int.out + * + * \note Even though the returned expression has dynamic size, in the case + * when it is applied to a fixed-size vector, it inherits a fixed maximal size, + * which means that evaluating it does not cause a dynamic memory allocation. + * + * \sa class Block, segment(int) + */ +template +inline VectorBlock MatrixBase + ::segment(int start, int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), start, size); +} + +/** This is the const version of segment(int,int).*/ +template +inline const VectorBlock +MatrixBase::segment(int start, int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), start, size); +} + +/** \returns a dynamic-size expression of the first coefficients of *this. + * + * \only_for_vectors + * + * \param size the number of coefficients in the block + * + * \addexample BlockInt \label How to reference a sub-vector (fixed-size) + * + * Example: \include MatrixBase_start_int.cpp + * Output: \verbinclude MatrixBase_start_int.out + * + * \note Even though the returned expression has dynamic size, in the case + * when it is applied to a fixed-size vector, it inherits a fixed maximal size, + * which means that evaluating it does not cause a dynamic memory allocation. + * + * \sa class Block, block(int,int) + */ +template +inline VectorBlock +MatrixBase::start(int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), 0, size); +} + +/** This is the const version of start(int).*/ +template +inline const VectorBlock +MatrixBase::start(int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), 0, size); +} + +/** \returns a dynamic-size expression of the last coefficients of *this. + * + * \only_for_vectors + * + * \param size the number of coefficients in the block + * + * \addexample BlockEnd \label How to reference the end of a vector (fixed-size) + * + * Example: \include MatrixBase_end_int.cpp + * Output: \verbinclude MatrixBase_end_int.out + * + * \note Even though the returned expression has dynamic size, in the case + * when it is applied to a fixed-size vector, it inherits a fixed maximal size, + * which means that evaluating it does not cause a dynamic memory allocation. + * + * \sa class Block, block(int,int) + */ +template +inline VectorBlock +MatrixBase::end(int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), this->size() - size, size); +} + +/** This is the const version of end(int).*/ +template +inline const VectorBlock +MatrixBase::end(int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), this->size() - size, size); +} + +/** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this + * + * \only_for_vectors + * + * The template parameter \a Size is the number of coefficients in the block + * + * \param start the index of the first element of the sub-vector + * + * Example: \include MatrixBase_template_int_segment.cpp + * Output: \verbinclude MatrixBase_template_int_segment.out + * + * \sa class Block + */ +template +template +inline VectorBlock +MatrixBase::segment(int start) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), start); +} + +/** This is the const version of segment(int).*/ +template +template +inline const VectorBlock +MatrixBase::segment(int start) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), start); +} + +/** \returns a fixed-size expression of the first coefficients of *this. + * + * \only_for_vectors + * + * The template parameter \a Size is the number of coefficients in the block + * + * \addexample BlockStart \label How to reference the start of a vector (fixed-size) + * + * Example: \include MatrixBase_template_int_start.cpp + * Output: \verbinclude MatrixBase_template_int_start.out + * + * \sa class Block + */ +template +template +inline VectorBlock +MatrixBase::start() +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), 0); +} + +/** This is the const version of start().*/ +template +template +inline const VectorBlock +MatrixBase::start() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), 0); +} + +/** \returns a fixed-size expression of the last coefficients of *this. + * + * \only_for_vectors + * + * The template parameter \a Size is the number of coefficients in the block + * + * Example: \include MatrixBase_template_int_end.cpp + * Output: \verbinclude MatrixBase_template_int_end.out + * + * \sa class Block + */ +template +template +inline VectorBlock +MatrixBase::end() +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), size() - Size); +} + +/** This is the const version of end.*/ +template +template +inline const VectorBlock +MatrixBase::end() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), size() - Size); +} + + +#endif // EIGEN_VECTORBLOCK_H diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 9606771a1..b457268af 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -39,6 +39,7 @@ template class Minor; template::Flags&DirectAccessBit ? DirectAccessBit : ei_traits::Flags&SparseBit> class Block; +template class VectorBlock; template class Transpose; template class Conjugate; template class CwiseNullaryOp; diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 7b139b0c1..b8d6aeb35 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -237,8 +237,6 @@ struct ei_special_scalar_op_base * TODO: could be a good idea to define a big ReturnType struct ?? */ template struct BlockReturnType { - typedef Block::RowsAtCompileTime == 1 ? 1 : RowsOrSize), - (ei_traits::ColsAtCompileTime == 1 ? 1 : RowsOrSize)> SubVectorType; typedef Block Type; }; diff --git a/doc/examples/class_FixedVectorBlock.cpp b/doc/examples/class_FixedVectorBlock.cpp new file mode 100644 index 000000000..c176be495 --- /dev/null +++ b/doc/examples/class_FixedVectorBlock.cpp @@ -0,0 +1,26 @@ +#include +USING_PART_OF_NAMESPACE_EIGEN +using namespace std; + +template +Eigen::VectorBlock +firstTwo(MatrixBase& v) +{ + return Eigen::VectorBlock(v.derived(), 0); +} + +template +const Eigen::VectorBlock +firstTwo(const MatrixBase& v) +{ + return Eigen::VectorBlock(v.derived(), 0); +} + +int main(int, char**) +{ + Matrix v; v << 1,2,3,4,5,6; + cout << firstTwo(4*v) << endl; // calls the const version + firstTwo(v) *= 2; // calls the non-const version + cout << "Now the vector v is:" << endl << v << endl; + return 0; +} diff --git a/doc/examples/class_VectorBlock.cpp b/doc/examples/class_VectorBlock.cpp new file mode 100644 index 000000000..d979f973a --- /dev/null +++ b/doc/examples/class_VectorBlock.cpp @@ -0,0 +1,26 @@ +#include +USING_PART_OF_NAMESPACE_EIGEN +using namespace std; + +template +Eigen::VectorBlock +segmentFromRange(MatrixBase& v, int start, int end) +{ + return Eigen::VectorBlock(v.derived(), start, end-start); +} + +template +const Eigen::VectorBlock +segmentFromRange(const MatrixBase& v, int start, int end) +{ + return Eigen::VectorBlock(v.derived(), start, end-start); +} + +int main(int, char**) +{ + Matrix v; v << 1,2,3,4,5,6; + cout << segmentFromRange(2*v, 2, 4) << endl; // calls the const version + segmentFromRange(v, 1, 3) *= 5; // calls the non-const version + cout << "Now the vector v is:" << endl << v << endl; + return 0; +} From ecc4f07af51f84959b8930c633f6e447a1d316f5 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 6 Jul 2009 09:17:25 +0200 Subject: [PATCH 11/25] fix doc of Quaternion::setFromTwoVectors --- Eigen/src/Geometry/Quaternion.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index d7480bf9c..f2ebece9a 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -342,11 +342,15 @@ Quaternion::toRotationMatrix(void) const return res; } -/** Sets *this to be a quaternion representing a rotation sending the vector \a a to the vector \a b. +/** Sets \c *this to be a quaternion representing a rotation between + * the two arbitrary vectors \a a and \a b. In other words, the built + * rotation represent a rotation sending the line of direction \a a + * to the line of direction \a b, both lines passing through the origin. * - * \returns a reference to *this. + * \returns a reference to \c *this. * - * Note that the two input vectors do \b not have to be normalized. + * Note that the two input vectors do \b not have to be normalized, and + * do not need to have the same norm. */ template template From 90f1e24579b8fa0605f63b2c9bfec5d023ca23ee Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 6 Jul 2009 10:35:20 +0200 Subject: [PATCH 12/25] significantly improve the accuracy of setFromTwoVectors (fixes #21) --- Eigen/Geometry | 1 + Eigen/src/Geometry/Quaternion.h | 24 ++++++++++++++---------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/Eigen/Geometry b/Eigen/Geometry index 7860d8fca..8698dd1a5 100644 --- a/Eigen/Geometry +++ b/Eigen/Geometry @@ -6,6 +6,7 @@ #include "src/Core/util/DisableMSVCWarnings.h" #include "Array" +#include "QR" #include #ifndef M_PI diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index f2ebece9a..088dff9f2 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -360,18 +360,22 @@ inline Quaternion& Quaternion::setFromTwoVectors(const MatrixBas Vector3 v1 = b.normalized(); Scalar c = v0.dot(v1); - // if dot == 1, vectors are the same - if (ei_isApprox(c,Scalar(1))) - { - // set to identity - this->w() = 1; this->vec().setZero(); - return *this; - } - // if dot == -1, vectors are opposites + // if dot == -1, vectors are nearly opposites + // => accuraletly compute the rotation axis by computing the + // intersection of the two planes. This is done by solving: + // x^T v0 = 0 + // x^T v1 = 0 + // under the constraint: + // ||x|| = 1 + // which yields an eigenvalue problem (or a SVD) if (ei_isApprox(c,Scalar(-1))) { - this->vec() = v0.unitOrthogonal(); - this->w() = 0; + c = std::max(c,-1); + SelfAdjointEigenSolver > eig(v0 * v0.transpose() + v1 * v1.transpose()); + Vector3 axis = eig.eigenvectors().col(0); + Scalar w2 = (Scalar(1)+c)*Scalar(0.5); + this->w() = ei_sqrt(w2); + this->vec() = axis * ei_sqrt(Scalar(1) - w2); return *this; } From 0cd158820cb8acb18507158fc1e3be327cdd0213 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 6 Jul 2009 11:15:38 +0200 Subject: [PATCH 13/25] switch from eigensolver to SVD which seems to be more accurate with float --- Eigen/Geometry | 2 +- Eigen/src/Geometry/Quaternion.h | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/Eigen/Geometry b/Eigen/Geometry index 8698dd1a5..9cae3459c 100644 --- a/Eigen/Geometry +++ b/Eigen/Geometry @@ -6,7 +6,7 @@ #include "src/Core/util/DisableMSVCWarnings.h" #include "Array" -#include "QR" +#include "SVD" #include #ifndef M_PI diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 088dff9f2..a76ccbdaf 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -367,12 +367,14 @@ inline Quaternion& Quaternion::setFromTwoVectors(const MatrixBas // x^T v1 = 0 // under the constraint: // ||x|| = 1 - // which yields an eigenvalue problem (or a SVD) + // which yields a singular value problem if (ei_isApprox(c,Scalar(-1))) { c = std::max(c,-1); - SelfAdjointEigenSolver > eig(v0 * v0.transpose() + v1 * v1.transpose()); - Vector3 axis = eig.eigenvectors().col(0); + + SVD > svd(v0 * v0.transpose() + v1 * v1.transpose()); + Vector3 axis = svd.matrixV().col(2); + Scalar w2 = (Scalar(1)+c)*Scalar(0.5); this->w() = ei_sqrt(w2); this->vec() = axis * ei_sqrt(Scalar(1) - w2); From 0c2232e5d972986ed90c917b68fb24eef372841b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 6 Jul 2009 13:47:41 +0200 Subject: [PATCH 14/25] quick reimplementation of SVD from the numeral recipes book: this is still not Eigen style code but at least it works for n>m and it is more accurate than the JAMA based version. (I needed it now, this is why I did that) --- Eigen/src/Geometry/Quaternion.h | 5 +- Eigen/src/SVD/SVD.h | 597 +++++++++++++------------------- test/svd.cpp | 2 +- 3 files changed, 254 insertions(+), 350 deletions(-) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index a76ccbdaf..8d1bbf9d2 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -371,13 +371,14 @@ inline Quaternion& Quaternion::setFromTwoVectors(const MatrixBas if (ei_isApprox(c,Scalar(-1))) { c = std::max(c,-1); - - SVD > svd(v0 * v0.transpose() + v1 * v1.transpose()); + Matrix m; m << v0.transpose(), v1.transpose(); + SVD > svd(m); Vector3 axis = svd.matrixV().col(2); Scalar w2 = (Scalar(1)+c)*Scalar(0.5); this->w() = ei_sqrt(w2); this->vec() = axis * ei_sqrt(Scalar(1) - w2); + return *this; } diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 45a7fbfa5..71f6763a8 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -34,9 +34,7 @@ * * \param MatrixType the type of the matrix of which we are computing the SVD decomposition * - * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N - * with \c M \>= \c N. - * + * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N. * * \sa MatrixBase::SVD() */ @@ -55,13 +53,13 @@ template class SVD typedef Matrix ColVector; typedef Matrix RowVector; - typedef Matrix MatrixUType; + typedef Matrix MatrixUType; typedef Matrix MatrixVType; - typedef Matrix SingularValuesType; + typedef Matrix SingularValuesType; public: - /** + /** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to @@ -70,9 +68,9 @@ template class SVD SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {} SVD(const MatrixType& matrix) - : m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())), + : m_matU(matrix.rows(), matrix.rows()), m_matV(matrix.cols(),matrix.cols()), - m_sigma(std::min(matrix.rows(),matrix.cols())), + m_sigma(matrix.cols()), m_isInitialized(false) { compute(matrix); @@ -81,22 +79,22 @@ template class SVD template bool solve(const MatrixBase &b, ResultType* result) const; - const MatrixUType& matrixU() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matU; - } - - const SingularValuesType& singularValues() const + const MatrixUType& matrixU() const { ei_assert(m_isInitialized && "SVD is not initialized."); - return m_sigma; + return m_matU; } - const MatrixVType& matrixV() const + const SingularValuesType& singularValues() const { ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matV; + return m_sigma; + } + + const MatrixVType& matrixV() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_matV; } void compute(const MatrixType& matrix); @@ -111,6 +109,23 @@ template class SVD template void computeScalingRotation(ScalingType *positive, RotationType *unitary) const; + protected: + // Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow. + inline static Scalar pythagora(Scalar a, Scalar b) + { + Scalar abs_a = ei_abs(a); + Scalar abs_b = ei_abs(b); + if (abs_a > abs_b) + return abs_a*ei_sqrt(Scalar(1.0)+ei_abs2(abs_b/abs_a)); + else + return (abs_b == Scalar(0.0) ? Scalar(0.0) : abs_b*ei_sqrt(Scalar(1.0)+ei_abs2(abs_a/abs_b))); + } + + inline static Scalar sign(Scalar a, Scalar b) + { + return (b >= Scalar(0.0) ? ei_abs(a) : -ei_abs(a)); + } + protected: /** \internal */ MatrixUType m_matU; @@ -123,7 +138,7 @@ template class SVD /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix * - * \note this code has been adapted from JAMA (public domain) + * \note this code has been adapted from Numerical Recipes, second edition. */ template void SVD::compute(const MatrixType& matrix) @@ -132,371 +147,259 @@ void SVD::compute(const MatrixType& matrix) const int n = matrix.cols(); const int nu = std::min(m,n); - m_matU.resize(m, nu); + m_matU.resize(m, m); m_matU.setZero(); - m_sigma.resize(std::min(m,n)); + m_sigma.resize(n); m_matV.resize(n,n); - RowVector e(n); - ColVector work(m); - MatrixType matA(matrix); - const bool wantu = true; - const bool wantv = true; - int i=0, j=0, k=0; + int max_iters = 30; - // Reduce A to bidiagonal form, storing the diagonal elements - // in s and the super-diagonal elements in e. - int nct = std::min(m-1,n); - int nrt = std::max(0,std::min(n-2,m)); - for (k = 0; k < std::max(nct,nrt); ++k) + MatrixVType& V = m_matV; + MatrixType A = matrix; + SingularValuesType& W = m_sigma; + + int flag,i,its,j,jj,k,l,nm; + Scalar anorm, c, f, g, h, s, scale, x, y, z; + bool convergence = true; + + Matrix rv1(n); + g = scale = anorm = 0; + // Householder reduction to bidiagonal form. + for (i=0; i= 0; k--) - { - if (m_sigma[k] != 0.0) - { - for (j = k+1; j < nu; ++j) + for (k=i; k0) - m_matU.col(k).start(k-1).setZero(); - } - else - { - m_matU.col(k).setZero(); - m_matU(k,k) = 1.0; + f = A(i, i); + g = -sign( ei_sqrt(s), f ); + h = f*g - s; + A(i, i)=f-g; + for (j=l; j= 0; k--) + W[i] = scale *g; + g = s = scale = 0.0; + if (i < m && i != (n-1)) { - if ((k < nrt) & (e[k] != 0.0)) + scale = A.row(i).end(n-l).cwise().abs().sum(); + if (scale) { - for (j = k+1; j < nu; ++j) + for (k=l; k=0; i--) + { + //Accumulation of right-hand transformations. + if (i < (n-1)) + { + if (g) + { + for (j=l; j::ret ? Scalar(-23) : Scalar(-52)); - while (p > 0) + // Accumulation of left-hand transformations. + for (i=std::min(m,n)-1; i>=0; i--) { - int k=0; - int kase=0; - - // Here is where a test for too many iterations would go. - - // This section of the program inspects for - // negligible elements in the s and e arrays. On - // completion the variables kase and k are set as follows. - - // kase = 1 if s(p) and e[k-1] are negligible and k

= -1; --k) + l = i+1; + g = W[i]; + for (j=l; j=0; k--) + { + for (its=1; its<=max_iters; its++) { - int ks; - for (ks = p-1; ks >= k; --ks) + flag=1; + for (l=k; l>=0; l--) { - if (ks == k) - break; - Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0)); - if (ei_abs(m_sigma[ks]) <= eps*t) + // Test for splitting. + nm=l-1; + // Note that rv1[1] is always zero. + if ((double)(ei_abs(rv1[l])+anorm) == anorm) { - m_sigma[ks] = 0.0; + flag=0; break; } + if ((double)(ei_abs(W[nm])+anorm) == anorm) + break; } - if (ks == k) + if (flag) { - kase = 3; + c=0.0; //Cancellation of rv1[l], if l > 1. + s=1.0; + for (i=l ;i<=k; i++) + { + f = s*rv1[i]; + rv1[i] = c*rv1[i]; + if ((double)(ei_abs(f)+anorm) == anorm) + break; + g = W[i]; + h = pythagora(f,g); + W[i] = h; + h = (Scalar)1.0/h; + c = g*h; + s = -f*h; + for (j=0; j= k; --j) - { - Scalar t(ei_hypot(m_sigma[j],f)); - Scalar cs(m_sigma[j]/t); - Scalar sn(f/t); - m_sigma[j] = t; - if (j != k) - { - f = -sn*e[j-1]; - e[j-1] = cs*e[j-1]; - } - if (wantv) - { - for (i = 0; i < n; ++i) - { - t = cs*m_matV(i,j) + sn*m_matV(i,p-1); - m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1); - m_matV(i,j) = t; - } - } - } - } - break; - - // Split at negligible s(k). - case 2: - { - Scalar f(e[k-1]); - e[k-1] = 0.0; - for (j = k; j < p; ++j) - { - Scalar t(ei_hypot(m_sigma[j],f)); - Scalar cs( m_sigma[j]/t); - Scalar sn(f/t); - m_sigma[j] = t; - f = -sn*e[j]; - e[j] = cs*e[j]; - if (wantu) - { - for (i = 0; i < m; ++i) - { - t = cs*m_matU(i,j) + sn*m_matU(i,k-1); - m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1); - m_matU(i,j) = t; - } - } - } - } - break; - - // Perform one qr step. - case 3: - { - // Calculate the shift. - Scalar scale = std::max(std::max(std::max(std::max( - ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])), - ei_abs(m_sigma[k])),ei_abs(e[k])); - Scalar sp = m_sigma[p-1]/scale; - Scalar spm1 = m_sigma[p-2]/scale; - Scalar epm1 = e[p-2]/scale; - Scalar sk = m_sigma[k]/scale; - Scalar ek = e[k]/scale; - Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2); - Scalar c = (sp*epm1)*(sp*epm1); - Scalar shift = 0.0; - if ((b != 0.0) || (c != 0.0)) - { - shift = ei_sqrt(b*b + c); - if (b < 0.0) - shift = -shift; - shift = c/(b + shift); - } - Scalar f = (sk + sp)*(sk - sp) + shift; - Scalar g = sk*ek; - - // Chase zeros. - - for (j = k; j < p-1; ++j) - { - Scalar t = ei_hypot(f,g); - Scalar cs = f/t; - Scalar sn = g/t; - if (j != k) - e[j-1] = t; - f = cs*m_sigma[j] + sn*e[j]; - e[j] = cs*e[j] - sn*m_sigma[j]; - g = sn*m_sigma[j+1]; - m_sigma[j+1] = cs*m_sigma[j+1]; - if (wantv) - { - for (i = 0; i < n; ++i) - { - t = cs*m_matV(i,j) + sn*m_matV(i,j+1); - m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1); - m_matV(i,j) = t; - } - } - t = ei_hypot(f,g); - cs = f/t; - sn = g/t; - m_sigma[j] = t; - f = cs*e[j] + sn*m_sigma[j+1]; - m_sigma[j+1] = -sn*e[j] + cs*m_sigma[j+1]; - g = sn*e[j+1]; - e[j+1] = cs*e[j+1]; - if (wantu && (j < m-1)) - { - for (i = 0; i < m; ++i) - { - t = cs*m_matU(i,j) + sn*m_matU(i,j+1); - m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1); - m_matU(i,j) = t; - } - } - } - e[p-2] = f; - iter = iter + 1; - } - break; - - // Convergence. - case 4: - { - // Make the singular values positive. - if (m_sigma[k] <= 0.0) - { - m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0); - if (wantv) - m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1); - } - - // Order the singular values. - while (k < pp) - { - if (m_sigma[k] >= m_sigma[k+1]) - break; - Scalar t = m_sigma[k]; - m_sigma[k] = m_sigma[k+1]; - m_sigma[k+1] = t; - if (wantv && (k < n-1)) - m_matV.col(k).swap(m_matV.col(k+1)); - if (wantu && (k < m-1)) - m_matU.col(k).swap(m_matU.col(k+1)); - ++k; - } - iter = 0; - p--; - } - break; - } // end big switch - } // end iterations + } + m_matU.setZero(); + if (m>=n) + m_matU.block(0,0,m,n) = A; + else + m_matU = A.block(0,0,m,m); m_isInitialized = true; } diff --git a/test/svd.cpp b/test/svd.cpp index 2eb7881a0..93e5ea017 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -50,7 +50,7 @@ template void svd(const MatrixType& m) MatrixType sigma = MatrixType::Zero(rows,cols); MatrixType matU = MatrixType::Zero(rows,rows); sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal(); - matU.block(0,0,rows,cols) = svd.matrixU(); + matU = svd.matrixU(); VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); } From f84bd3e7b1b669c9ad1469dbdd2929531db6645f Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 6 Jul 2009 15:01:30 +0200 Subject: [PATCH 15/25] include the fixes of the third edition --- Eigen/src/SVD/SVD.h | 129 ++++++++++++++++++++++---------------------- 1 file changed, 66 insertions(+), 63 deletions(-) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 71f6763a8..f27b7e66d 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -111,7 +111,7 @@ template class SVD protected: // Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow. - inline static Scalar pythagora(Scalar a, Scalar b) + inline static Scalar pythag(Scalar a, Scalar b) { Scalar abs_a = ei_abs(a); Scalar abs_b = ei_abs(b); @@ -138,14 +138,13 @@ template class SVD /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix * - * \note this code has been adapted from Numerical Recipes, second edition. + * \note this code has been adapted from Numerical Recipes, third edition. */ template void SVD::compute(const MatrixType& matrix) { const int m = matrix.rows(); const int n = matrix.cols(); - const int nu = std::min(m,n); m_matU.resize(m, m); m_matU.setZero(); @@ -158,22 +157,24 @@ void SVD::compute(const MatrixType& matrix) MatrixType A = matrix; SingularValuesType& W = m_sigma; - int flag,i,its,j,jj,k,l,nm; + bool flag; + int i,its,j,jj,k,l,nm; Scalar anorm, c, f, g, h, s, scale, x, y, z; bool convergence = true; + Scalar eps = precision(); Matrix rv1(n); g = scale = anorm = 0; // Householder reduction to bidiagonal form. for (i=0; i::compute(const MatrixType& matrix) g = -sign( ei_sqrt(s), f ); h = f*g - s; A(i, i)=f-g; - for (j=l; j::compute(const MatrixType& matrix) A.col(i).end(m-i) *= scale; } } - W[i] = scale *g; + W[i] = scale * g; g = s = scale = 0.0; - if (i < m && i != (n-1)) + if (i+1 <= m && i+1 != n) { - scale = A.row(i).end(n-l).cwise().abs().sum(); - if (scale) + scale = A.row(i).end(n-l+1).cwise().abs().sum(); + if (scale != Scalar(0)) { - for (k=l; k=0; i--) + for (i=n-1; i>=0; i--) { //Accumulation of right-hand transformations. - if (i < (n-1)) + if (i < n-1) { - if (g) + if (g != Scalar(0.0)) { - for (j=l; j::compute(const MatrixType& matrix) { l = i+1; g = W[i]; - for (j=l; j0) + A.row(i).end(n-l).setZero(); + if (g != Scalar(0.0)) { - g = (Scalar)1.0/g; + g = Scalar(1.0)/g; for (j=l; j=0; k--) + for (k=n-1; k>=0; k--) { - for (its=1; its<=max_iters; its++) + for (its=0; its=0; l--) { // Test for splitting. - nm=l-1; + nm = l-1; // Note that rv1[1] is always zero. - if ((double)(ei_abs(rv1[l])+anorm) == anorm) + //if ((double)(ei_abs(rv1[l])+anorm) == anorm) + if (l==0 || ei_abs(rv1[l]) <= eps*anorm) { - flag=0; + flag = false; break; } - if ((double)(ei_abs(W[nm])+anorm) == anorm) + //if ((double)(ei_abs(W[nm])+anorm) == anorm) + if (ei_abs(W[nm]) <= eps*anorm) break; } if (flag) { - c=0.0; //Cancellation of rv1[l], if l > 1. - s=1.0; - for (i=l ;i<=k; i++) + c = 0.0; //Cancellation of rv1[l], if l > 0. + s = 1.0; + for (i=l ;i::compute(const MatrixType& matrix) } break; } - if (its == max_iters) + if (its+1 == max_iters) { convergence = false; } @@ -329,19 +332,19 @@ void SVD::compute(const MatrixType& matrix) y = W[nm]; g = rv1[nm]; h = rv1[k]; - f = ((y-z)*(y+z) + (g-h)*(g+h))/((Scalar)2.0*h*y); - g = pythagora(f,1.0); + f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y); + g = pythag(f,1.0); f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x; c = s = 1.0; //Next QR transformation: - for (j=l; j<= nm;j++) + for (j=l; j<=nm; j++) { i = j+1; g = rv1[i]; y = W[i]; h = s*g; g = c*g; - z = pythagora(f,h); + z = pythag(f,h); rv1[j] = z; c = f/z; s = h/z; @@ -351,15 +354,15 @@ void SVD::compute(const MatrixType& matrix) y *= c; for (jj=0; jj::compute(const MatrixType& matrix) x = c*y - s*g; for (jj=0; jj Date: Mon, 6 Jul 2009 17:12:10 +0200 Subject: [PATCH 16/25] * rename QR to HouseholderQR because really that impacts the API, not just the impl. * rename qr() to householderQr(), for same reason. * clarify that it's non-pivoting, non-rank-revealing, so remove all the rank API, make solve() be void instead of bool, update the docs/test, etc. * fix warning in SVD --- Eigen/src/Core/MatrixBase.h | 2 +- Eigen/src/Core/util/ForwardDeclarations.h | 2 +- Eigen/src/QR/QR.h | 161 ++++-------------- Eigen/src/SVD/SVD.h | 1 - doc/eigendoxy_header.html.in | 2 +- .../{QR_solve.cpp => HouseholderQR_solve.cpp} | 11 +- test/geo_hyperplane.cpp | 2 +- test/main.h | 4 +- test/qr.cpp | 62 +------ 9 files changed, 43 insertions(+), 204 deletions(-) rename doc/snippets/{QR_solve.cpp => HouseholderQR_solve.cpp} (52%) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 65ab02d62..941539214 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -653,7 +653,7 @@ template class MatrixBase /////////// QR module /////////// - const QR qr() const; + const HouseholderQR householderQr() const; EigenvaluesReturnType eigenvalues() const; RealScalar operatorNorm() const; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index b457268af..a2105604a 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -112,7 +112,7 @@ template class Reverse; template class LU; template class PartialLU; -template class QR; +template class HouseholderQR; template class SVD; template class LLT; template class LDLT; diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 90751dd42..d6d3d2081 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -22,24 +22,26 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#ifndef EIGEN_QR_H -#define EIGEN_QR_H +#ifndef EIGEN_HouseholderQR_H +#define EIGEN_HouseholderQR_H -/** \ingroup QR_Module +/** \ingroup HouseholderQR_Module * \nonstableyet * - * \class QR + * \class HouseholderQR * - * \brief QR decomposition of a matrix + * \brief Householder QR decomposition of a matrix * * \param MatrixType the type of the matrix of which we are computing the QR decomposition * * This class performs a QR decomposition using Householder transformations. The result is * stored in a compact way compatible with LAPACK. * + * Note that no pivoting is performed. This is \b not a rank-revealing decomposition. + * * \sa MatrixBase::qr() */ -template class QR +template class HouseholderQR { public: @@ -53,88 +55,23 @@ template class QR * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to - * perform decompositions via QR::compute(const MatrixType&). + * perform decompositions via HouseholderQR::compute(const MatrixType&). */ - QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} - QR(const MatrixType& matrix) + HouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs(matrix.cols()), m_isInitialized(false) { compute(matrix); } - - /** \deprecated use isInjective() - * \returns whether or not the matrix is of full rank - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - EIGEN_DEPRECATED bool isFullRank() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.cols(); - } - - /** \returns the rank of the matrix of which *this is the QR decomposition. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - int rank() const; - - /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline int dimensionOfKernel() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return m_qr.cols() - rank(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents an injective - * linear map, i.e. has trivial kernel; false otherwise. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isInjective() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.cols(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents a surjective - * linear map; false otherwise. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isSurjective() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.rows(); - } - - /** \returns true if the matrix of which *this is the QR decomposition is invertible. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isInvertible() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return isInjective() && isSurjective(); - } - + /** \returns a read-only expression of the matrix R of the actual the QR decomposition */ const Part, UpperTriangular> matrixR(void) const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); int cols = m_qr.cols(); return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part(); } @@ -148,22 +85,14 @@ template class QR * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). * If no solution exists, *result is left with undefined coefficients. * - * \returns true if any solution exists, false if no solution exists. - * - * \note If there exist more than one solution, this method will arbitrarily choose one. - * If you need a complete analysis of the space of solutions, take the one solution obtained - * by this method and add to it elements of the kernel, as determined by kernel(). - * * \note The case where b is a matrix is not yet implemented. Also, this * code is space inefficient. * - * Example: \include QR_solve.cpp - * Output: \verbinclude QR_solve.out - * - * \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse() + * Example: \include HouseholderQR_solve.cpp + * Output: \verbinclude HouseholderQR_solve.out */ template - bool solve(const MatrixBase& b, ResultType *result) const; + void solve(const MatrixBase& b, ResultType *result) const; MatrixType matrixQ(void) const; @@ -172,34 +101,14 @@ template class QR protected: MatrixType m_qr; VectorType m_hCoeffs; - mutable int m_rank; - mutable bool m_rankIsUptodate; bool m_isInitialized; }; -/** \returns the rank of the matrix of which *this is the QR decomposition. */ -template -int QR::rank() const -{ - ei_assert(m_isInitialized && "QR is not initialized."); - if (!m_rankIsUptodate) - { - RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff(); - int n = m_qr.cols(); - m_rank = 0; - while(m_rank -void QR::compute(const MatrixType& matrix) +void HouseholderQR::compute(const MatrixType& matrix) { - m_rankIsUptodate = false; m_qr = matrix; m_hCoeffs.resize(matrix.cols()); @@ -262,12 +171,12 @@ void QR::compute(const MatrixType& matrix) template template -bool QR::solve( +void HouseholderQR::solve( const MatrixBase& b, ResultType *result ) const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); const int rows = m_qr.rows(); ei_assert(b.rows() == rows); result->resize(rows, b.cols()); @@ -276,27 +185,17 @@ bool QR::solve( // Q^T without explicitly forming matrixQ(). Investigate. *result = matrixQ().transpose()*b; - if(!isSurjective()) - { - // is result is in the image of R ? - RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff(); - for(int col = 0; col < result->cols(); ++col) - for(int row = m_rank; row < result->rows(); ++row) - if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res)) - return false; - } - m_qr.corner(TopLeft, m_rank, m_rank) + const int rank = std::min(result->rows(), result->cols()); + m_qr.corner(TopLeft, rank, rank) .template marked() - .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols())); - - return true; + .solveTriangularInPlace(result->corner(TopLeft, rank, result->cols())); } /** \returns the matrix Q */ template -MatrixType QR::matrixQ() const +MatrixType HouseholderQR::matrixQ() const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); // compute the product Q_0 Q_1 ... Q_n-1, // where Q_k is the k-th Householder transformation I - h_k v_k v_k' // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] @@ -319,16 +218,16 @@ MatrixType QR::matrixQ() const #endif // EIGEN_HIDE_HEAVY_CODE -/** \return the QR decomposition of \c *this. +/** \return the Householder QR decomposition of \c *this. * - * \sa class QR + * \sa class HouseholderQR */ template -const QR::PlainMatrixType> -MatrixBase::qr() const +const HouseholderQR::PlainMatrixType> +MatrixBase::householderQr() const { - return QR(eval()); + return HouseholderQR(eval()); } -#endif // EIGEN_QR_H +#endif // EIGEN_HouseholderQR_H diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 71f6763a8..97f82c78f 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -145,7 +145,6 @@ void SVD::compute(const MatrixType& matrix) { const int m = matrix.rows(); const int n = matrix.cols(); - const int nu = std::min(m,n); m_matU.resize(m, m); m_matU.setZero(); diff --git a/doc/eigendoxy_header.html.in b/doc/eigendoxy_header.html.in index 97cafd41d..572c47158 100644 --- a/doc/eigendoxy_header.html.in +++ b/doc/eigendoxy_header.html.in @@ -6,5 +6,5 @@ \ No newline at end of file diff --git a/doc/snippets/QR_solve.cpp b/doc/snippets/HouseholderQR_solve.cpp similarity index 52% rename from doc/snippets/QR_solve.cpp rename to doc/snippets/HouseholderQR_solve.cpp index 7e629f851..aa9404951 100644 --- a/doc/snippets/QR_solve.cpp +++ b/doc/snippets/HouseholderQR_solve.cpp @@ -4,11 +4,6 @@ Matrix3f y = Matrix3f::Random(); cout << "Here is the matrix m:" << endl << m << endl; cout << "Here is the matrix y:" << endl << y << endl; Matrix3f x; -if(m.qr().solve(y, &x)) -{ - assert(y.isApprox(m*x)); - cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; -} -else - cout << "The equation mx=y does not have any solution." << endl; - +m.householderQr().solve(y, &x)); +assert(y.isApprox(m*x)); +cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp index c5d2715db..f8281a16b 100644 --- a/test/geo_hyperplane.cpp +++ b/test/geo_hyperplane.cpp @@ -64,7 +64,7 @@ template void hyperplane(const HyperplaneType& _plane) // transform if (!NumTraits::IsComplex) { - MatrixType rot = MatrixType::Random(dim,dim).qr().matrixQ(); + MatrixType rot = MatrixType::Random(dim,dim).householderQr().matrixQ(); DiagonalMatrix scaling(VectorType::Random()); Translation translation(VectorType::Random()); diff --git a/test/main.h b/test/main.h index a6a780603..53c8245c6 100644 --- a/test/main.h +++ b/test/main.h @@ -241,8 +241,8 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri 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); - QR qra(a); - QR qrb(b); + HouseholderQR qra(a); + HouseholderQR qrb(b); m = (qra.matrixQ() * d * qrb.matrixQ()).lazy(); } diff --git a/test/qr.cpp b/test/qr.cpp index 6e96f1e97..88a447c4b 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -36,7 +36,7 @@ template void qr(const MatrixType& m) typedef Matrix VectorType; MatrixType a = MatrixType::Random(rows,cols); - QR qrOfA(a); + HouseholderQR qrOfA(a); VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR()); VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR()); @@ -55,40 +55,6 @@ template void qr(const MatrixType& m) VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); } -template void qr_non_invertible() -{ - /* this test covers the following files: QR.h */ - int rows = ei_random(20,200), cols = ei_random(20,rows), cols2 = ei_random(20,rows); - int rank = ei_random(1, std::min(rows, cols)-1); - - MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); - createRandomMatrixOfRank(rank, rows, cols, m1); - - QR lu(m1); -// typename LU::KernelResultType m1kernel = lu.kernel(); -// typename LU::ImageResultType m1image = lu.image(); - std::cerr << rows << "x" << cols << " " << rank << " " << lu.rank() << "\n"; - if (rank != lu.rank()) - std::cerr << lu.matrixR().diagonal().transpose() << "\n"; - VERIFY(rank == lu.rank()); - VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); - VERIFY(!lu.isInjective()); - VERIFY(!lu.isInvertible()); - VERIFY(lu.isSurjective() == (lu.rank() == rows)); -// VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); -// VERIFY(m1image.lu().rank() == rank); -// MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); -// sidebyside << m1, m1image; -// VERIFY(sidebyside.lu().rank() == rank); - m2 = MatrixType::Random(cols,cols2); - m3 = m1*m2; - m2 = MatrixType::Random(cols,cols2); - lu.solve(m3, &m2); - VERIFY_IS_APPROX(m3, m1*m2); - m3 = MatrixType::Random(rows,cols2); - VERIFY(!lu.solve(m3, &m2)); -} - template void qr_invertible() { /* this test covers the following files: QR.h */ @@ -105,33 +71,18 @@ template void qr_invertible() m1 += a * a.adjoint(); } - QR lu(m1); - VERIFY(0 == lu.dimensionOfKernel()); - VERIFY(size == lu.rank()); - VERIFY(lu.isInjective()); - VERIFY(lu.isSurjective()); - VERIFY(lu.isInvertible()); -// VERIFY(lu.image().lu().isInvertible()); + HouseholderQR qr(m1); m3 = MatrixType::Random(size,size); - lu.solve(m3, &m2); + qr.solve(m3, &m2); //std::cerr << m3 - m1*m2 << "\n\n"; VERIFY_IS_APPROX(m3, m1*m2); -// VERIFY_IS_APPROX(m2, lu.inverse()*m3); - m3 = MatrixType::Random(size,size); - VERIFY(lu.solve(m3, &m2)); } template void qr_verify_assert() { MatrixType tmp; - QR qr; - VERIFY_RAISES_ASSERT(qr.isFullRank()) - VERIFY_RAISES_ASSERT(qr.rank()) - VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) - VERIFY_RAISES_ASSERT(qr.isInjective()) - VERIFY_RAISES_ASSERT(qr.isSurjective()) - VERIFY_RAISES_ASSERT(qr.isInvertible()) + HouseholderQR qr; VERIFY_RAISES_ASSERT(qr.matrixR()) VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) VERIFY_RAISES_ASSERT(qr.matrixQ()) @@ -149,11 +100,6 @@ void test_qr() } for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( qr_non_invertible() ); - CALL_SUBTEST( qr_non_invertible() ); - // TODO fix issue with complex -// CALL_SUBTEST( qr_non_invertible() ); -// CALL_SUBTEST( qr_non_invertible() ); CALL_SUBTEST( qr_invertible() ); CALL_SUBTEST( qr_invertible() ); // TODO fix issue with complex From e057beee4ee682f3128d36c74f3aa57fc9ce148e Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 6 Jul 2009 17:20:07 +0200 Subject: [PATCH 17/25] fix some search-and-replace damage --- Eigen/src/QR/QR.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index d6d3d2081..a83ec5ec8 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -22,10 +22,10 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#ifndef EIGEN_HouseholderQR_H -#define EIGEN_HouseholderQR_H +#ifndef EIGEN_QR_H +#define EIGEN_QR_H -/** \ingroup HouseholderQR_Module +/** \ingroup QR_Module * \nonstableyet * * \class HouseholderQR @@ -230,4 +230,4 @@ MatrixBase::householderQr() const } -#endif // EIGEN_HouseholderQR_H +#endif // EIGEN_QR_H From 889726bf7cf82b30e4a7140f467a175bab1dca2d Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 6 Jul 2009 17:24:11 +0200 Subject: [PATCH 18/25] add matrixQR() method exposing the storage. that's where the householder thing impacts the API. --- Eigen/src/QR/QR.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index a83ec5ec8..c5e8f9097 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -95,6 +95,11 @@ template class HouseholderQR void solve(const MatrixBase& b, ResultType *result) const; MatrixType matrixQ(void) const; + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + * in a LAPACK-compatible way. + */ + const MatrixType& matrixQR() const { return m_qr; } void compute(const MatrixType& matrix); From 629e083d81a9203f4b2bc07a6cf0f8a61eef07c7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 10 Jul 2009 08:21:20 +0200 Subject: [PATCH 19/25] slight change in the comparison to -1 --- Eigen/src/Geometry/Quaternion.h | 6 ++---- Eigen/src/SVD/SVD.h | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 8d1bbf9d2..9385f259d 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -368,20 +368,18 @@ inline Quaternion& Quaternion::setFromTwoVectors(const MatrixBas // under the constraint: // ||x|| = 1 // which yields a singular value problem - if (ei_isApprox(c,Scalar(-1))) + if (c < Scalar(-1)+precision()) { c = std::max(c,-1); Matrix m; m << v0.transpose(), v1.transpose(); SVD > svd(m); Vector3 axis = svd.matrixV().col(2); - + Scalar w2 = (Scalar(1)+c)*Scalar(0.5); this->w() = ei_sqrt(w2); this->vec() = axis * ei_sqrt(Scalar(1) - w2); - return *this; } - Vector3 axis = v0.cross(v1); Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2)); Scalar invs = Scalar(1)/s; diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index f27b7e66d..f9f9feb89 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2008-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public From ec5c608aa3caea7f14aa03e4e13041ee9e2664de Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 10 Jul 2009 16:10:03 +0200 Subject: [PATCH 20/25] Set of fixes and workaround to make sun studio more happy. Still remains the problem of alignment and vectorization. --- Eigen/Core | 2 + Eigen/src/Array/Replicate.h | 17 -------- Eigen/src/Array/VectorwiseOp.h | 52 +++++++++++++---------- Eigen/src/Core/Block.h | 6 +++ Eigen/src/Core/CwiseUnaryView.h | 13 +++--- Eigen/src/Core/MapBase.h | 1 + Eigen/src/Core/Matrix.h | 3 +- Eigen/src/Core/MatrixBase.h | 8 +++- Eigen/src/Core/Product.h | 15 +++---- Eigen/src/Core/util/Constants.h | 4 +- Eigen/src/Core/util/ForwardDeclarations.h | 2 +- Eigen/src/Core/util/Macros.h | 16 ++++--- cmake/FindEigen2.cmake | 2 +- test/array.cpp | 2 +- test/basicstuff.cpp | 10 ++--- 15 files changed, 81 insertions(+), 72 deletions(-) diff --git a/Eigen/Core b/Eigen/Core index 0d3650c99..a7a2a768a 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -88,6 +88,8 @@ #include #include #include +// for min/max: +#include #if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(EIGEN_NO_EXCEPTIONS) #define EIGEN_EXCEPTIONS diff --git a/Eigen/src/Array/Replicate.h b/Eigen/src/Array/Replicate.h index df3afbbdb..02f9c0601 100644 --- a/Eigen/src/Array/Replicate.h +++ b/Eigen/src/Array/Replicate.h @@ -140,21 +140,4 @@ VectorwiseOp::replicate(int factor) const (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); } -/** \nonstableyet - * \return an expression of the replication of each column (or row) of \c *this - * - * Example: \include DirectionWise_replicate.cpp - * Output: \verbinclude DirectionWise_replicate.out - * - * \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate - */ -template -template -const Replicate -VectorwiseOp::replicate(int factor) const -{ - return Replicate - (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); -} - #endif // EIGEN_REPLICATE_H diff --git a/Eigen/src/Array/VectorwiseOp.h b/Eigen/src/Array/VectorwiseOp.h index 50302bba4..3684d7cd1 100644 --- a/Eigen/src/Array/VectorwiseOp.h +++ b/Eigen/src/Array/VectorwiseOp.h @@ -179,6 +179,11 @@ template class VectorwiseOp > Type; }; + enum { + IsVertical = (Direction==Vertical) ? 1 : 0, + IsHorizontal = (Direction==Horizontal) ? 1 : 0 + }; + protected: /** \internal @@ -222,9 +227,17 @@ template class VectorwiseOp /** \internal */ inline const ExpressionType& _expression() const { return m_matrix; } + /** \returns a row or column vector expression of \c *this reduxed by \a func + * + * The template parameter \a BinaryOp is the type of the functor + * of the custom redux operator. Note that func must be an associative operator. + * + * \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise() + */ template const typename ReduxReturnType::Type - redux(const BinaryOp& func = BinaryOp()) const; + redux(const BinaryOp& func = BinaryOp()) const + { return typename ReduxReturnType::Type(_expression(), func); } /** \returns a row (or column) vector expression of the smallest coefficient * of each column (or row) of the referenced expression. @@ -319,16 +332,26 @@ template class VectorwiseOp * * \sa MatrixBase::reverse() */ const Reverse reverse() const - { - return Reverse( _expression() ); - } + { return Reverse( _expression() ); } const Replicate replicate(int factor) const; - template - const Replicate - replicate(int factor = Factor) const; + /** \nonstableyet + * \return an expression of the replication of each column (or row) of \c *this + * + * Example: \include DirectionWise_replicate.cpp + * Output: \verbinclude DirectionWise_replicate.out + * + * \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate + */ + // NOTE implemented here because of sunstudio's compilation errors + template const Replicate + replicate(int factor = Factor) const + { + return Replicate + (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); + } /////////// Artithmetic operators /////////// @@ -466,19 +489,4 @@ MatrixBase::rowwise() return derived(); } -/** \returns a row or column vector expression of \c *this reduxed by \a func - * - * The template parameter \a BinaryOp is the type of the functor - * of the custom redux operator. Note that func must be an associative operator. - * - * \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise() - */ -template -template -const typename VectorwiseOp::template ReduxReturnType::Type -VectorwiseOp::redux(const BinaryOp& func) const -{ - return typename ReduxReturnType::Type(_expression(), func); -} - #endif // EIGEN_PARTIAL_REDUX_H diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h index d7af971f9..382518696 100644 --- a/Eigen/src/Core/Block.h +++ b/Eigen/src/Core/Block.h @@ -271,13 +271,19 @@ class Block inline int stride(void) const { return m_matrix.stride(); } + #ifndef __SUNPRO_CC + // FIXME sunstudio is not friendly with the above friend... protected: + #endif + #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal used by allowAligned() */ inline Block(const MatrixType& matrix, const Scalar* data, int blockRows, int blockCols) : Base(data, blockRows, blockCols), m_matrix(matrix) {} + #endif + protected: const typename MatrixType::Nested m_matrix; }; diff --git a/Eigen/src/Core/CwiseUnaryView.h b/Eigen/src/Core/CwiseUnaryView.h index c8fd98ea1..8d8b51b4c 100644 --- a/Eigen/src/Core/CwiseUnaryView.h +++ b/Eigen/src/Core/CwiseUnaryView.h @@ -42,13 +42,13 @@ struct ei_traits > : ei_traits { typedef typename ei_result_of< - ViewOp(typename MatrixType::Scalar) + ViewOp(typename ei_traits::Scalar) >::type Scalar; typedef typename MatrixType::Nested MatrixTypeNested; typedef typename ei_unref::type _MatrixTypeNested; enum { - Flags = (_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), - CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits::Cost + Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), + CoeffReadCost = ei_traits<_MatrixTypeNested>::CoeffReadCost + ei_functor_traits::Cost }; }; @@ -62,7 +62,7 @@ class CwiseUnaryView : ei_no_assignment_operator, inline CwiseUnaryView(const MatrixType& mat, const ViewOp& func = ViewOp()) : m_matrix(mat), m_functor(func) {} - + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(CwiseUnaryView) EIGEN_STRONG_INLINE int rows() const { return m_matrix.rows(); } @@ -77,7 +77,7 @@ class CwiseUnaryView : ei_no_assignment_operator, { return m_functor(m_matrix.coeff(index)); } - + EIGEN_STRONG_INLINE Scalar& coeffRef(int row, int col) { return m_functor(m_matrix.const_cast_derived().coeffRef(row, col)); @@ -89,7 +89,8 @@ class CwiseUnaryView : ei_no_assignment_operator, } protected: - const typename MatrixType::Nested m_matrix; + // FIXME changed from MatrixType::Nested because of a weird compilation error with sun CC + const typename ei_nested::type m_matrix; const ViewOp m_functor; }; diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index 126974918..59bf69ad6 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -178,6 +178,7 @@ template class MapBase } using Base::operator*=; + using Base::operator+=; template Derived& operator+=(const MatrixBase& other) diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index e099ba1e7..5301f4849 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -124,6 +124,7 @@ class Matrix { public: EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix) + enum { Options = _Options }; friend class Eigen::Map; typedef class Eigen::Map UnalignedMapType; @@ -217,7 +218,7 @@ class Matrix * * This method is intended for dynamic-size matrices, although it is legal to call it on any * matrix as long as fixed dimensions are left unchanged. If you only want to change the number - * of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t). + * of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t). * * If the current number of coefficients of \c *this exactly matches the * product \a rows * \a cols, then no memory allocation is performed and diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 941539214..2a49ae620 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -137,10 +137,14 @@ template class MatrixBase * constructed from this one. See the \ref flags "list of flags". */ - CoeffReadCost = ei_traits::CoeffReadCost + CoeffReadCost = ei_traits::CoeffReadCost, /**< This is a rough measure of how expensive it is to read one coefficient from * this expression. */ + +#ifndef EIGEN_PARSED_BY_DOXYGEN + _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC +#endif }; /** Default constructor. Just checks at compile-time for self-consistency of the flags. */ @@ -204,7 +208,7 @@ template class MatrixBase /** \internal the return type of coeff() */ - typedef typename ei_meta_if::ret CoeffReturnType; + typedef typename ei_meta_if<_HasDirectAccess, const Scalar&, Scalar>::ret CoeffReturnType; /** \internal Represents a matrix with all coefficients equal to one another*/ typedef CwiseNullaryOp,Derived> ConstantReturnType; diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 353aec48e..2589d548e 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -83,15 +83,14 @@ struct ProductReturnType template struct ei_product_mode { enum{ - value = ei_is_diagonal::ret || ei_is_diagonal::ret ? DiagonalProduct - : Lhs::MaxColsAtCompileTime == Dynamic - && ( Lhs::MaxRowsAtCompileTime == Dynamic - || Rhs::MaxColsAtCompileTime == Dynamic ) - && (!(Rhs::IsVectorAtCompileTime && (Lhs::Flags&RowMajorBit) && (!(Lhs::Flags&DirectAccessBit)))) - && (!(Lhs::IsVectorAtCompileTime && (!(Rhs::Flags&RowMajorBit)) && (!(Rhs::Flags&DirectAccessBit)))) - && (ei_is_same_type::ret) + : ei_traits::MaxColsAtCompileTime == Dynamic + && ( ei_traits::MaxRowsAtCompileTime == Dynamic + || ei_traits::MaxColsAtCompileTime == Dynamic ) + && (!(ei_traits::IsVectorAtCompileTime && (ei_traits::Flags&RowMajorBit) && (!(ei_traits::Flags&DirectAccessBit)))) + && (!(ei_traits::IsVectorAtCompileTime && (!(ei_traits::Flags&RowMajorBit)) && (!(ei_traits::Flags&DirectAccessBit)))) + && (ei_is_same_type::Scalar, typename ei_traits::Scalar>::ret) ? CacheFriendlyProduct : NormalProduct }; }; @@ -215,7 +214,7 @@ template class Product */ EIGEN_STRONG_INLINE bool _useCacheFriendlyProduct() const { - return m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + return m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD && ( rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD || cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD); } diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index c2641537a..572cf7caf 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -188,7 +188,7 @@ const unsigned int HereditaryBits = RowMajorBit // diagonal means both upper and lower triangular const unsigned DiagonalBits = UpperTriangularBit | LowerTriangularBit; - + // Possible values for the Mode parameter of part() const unsigned int UpperTriangular = UpperTriangularBit; const unsigned int StrictlyUpperTriangular = UpperTriangularBit | ZeroDiagBit; @@ -201,7 +201,7 @@ const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit; template struct ei_is_diagonal { enum { - ret = ( (unsigned int)(T::Flags) & DiagonalBits ) == DiagonalBits + ret = ( int(ei_traits::Flags) & DiagonalBits ) == DiagonalBits }; }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index a2105604a..ad7256635 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -97,7 +97,7 @@ template struct ei_scalar_multiple2_op; struct IOFormat; template -void ei_cache_friendly_product( +static void ei_cache_friendly_product( int _rows, int _cols, int depth, bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride, bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride, diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 8f62f9b69..2b2444ef0 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -51,7 +51,8 @@ #define EIGEN_GCC3_OR_OLDER 0 #endif -#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER +// FIXME vectorization + alignment is completely disabled with sun studio +#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC) #define EIGEN_ARCH_WANTS_ALIGNMENT 1 #else #define EIGEN_ARCH_WANTS_ALIGNMENT 0 @@ -97,7 +98,7 @@ /** Allows to disable some optimizations which might affect the accuracy of the result. * Such optimization are enabled by default, and set EIGEN_FAST_MATH to 0 to disable them. * They currently include: - * - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled. + * - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled. */ #ifndef EIGEN_FAST_MATH #define EIGEN_FAST_MATH 1 @@ -199,13 +200,16 @@ using Eigen::ei_cos; * vectorized and non-vectorized code. */ #if !EIGEN_ALIGN -#define EIGEN_ALIGN_128 + #define EIGEN_ALIGN_128 #elif (defined __GNUC__) -#define EIGEN_ALIGN_128 __attribute__((aligned(16))) + #define EIGEN_ALIGN_128 __attribute__((aligned(16))) #elif (defined _MSC_VER) -#define EIGEN_ALIGN_128 __declspec(align(16)) + #define EIGEN_ALIGN_128 __declspec(align(16)) +#elif (defined __SUNPRO_CC) + // FIXME not sure about this one: + #define EIGEN_ALIGN_128 __attribute__((aligned(16))) #else -#error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler + #error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler #endif #define EIGEN_RESTRICT __restrict diff --git a/cmake/FindEigen2.cmake b/cmake/FindEigen2.cmake index 0e86aaf29..ee054b0da 100644 --- a/cmake/FindEigen2.cmake +++ b/cmake/FindEigen2.cmake @@ -46,7 +46,7 @@ macro(_eigen2_check_version) endif(${EIGEN2_VERSION} VERSION_LESS ${Eigen2_FIND_VERSION}) if(NOT EIGEN2_VERSION_OK) - + message(STATUS "Eigen2 version ${EIGEN2_VERSION} found in ${EIGEN2_INCLUDE_DIR}, " "but at least version ${Eigen2_FIND_VERSION} is required") endif(NOT EIGEN2_VERSION_OK) diff --git a/test/array.cpp b/test/array.cpp index 37deeaa4f..a308f7366 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2008-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 2f088988f..52aff93ee 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -107,7 +107,7 @@ template void basicStuff(const MatrixType& m) { VERIFY_IS_NOT_APPROX(m3, m1); } - + m3.real() = m1.real(); VERIFY_IS_APPROX(static_cast(m3).real(), static_cast(m1).real()); VERIFY_IS_APPROX(static_cast(m3).real(), m1.real()); @@ -121,16 +121,16 @@ template void basicStuffComplex(const MatrixType& m) int rows = m.rows(); int cols = m.cols(); - + Scalar s1 = ei_random(), s2 = ei_random(); - + VERIFY(ei_real(s1)==ei_real_ref(s1)); VERIFY(ei_imag(s1)==ei_imag_ref(s1)); ei_real_ref(s1) = ei_real(s2); ei_imag_ref(s1) = ei_imag(s2); VERIFY(s1==s2); - + RealMatrixType rm1 = RealMatrixType::Random(rows,cols), rm2 = RealMatrixType::Random(rows,cols); MatrixType cm(rows,cols); @@ -162,7 +162,7 @@ void test_basicstuff() CALL_SUBTEST( basicStuff(MatrixXcd(20, 20)) ); CALL_SUBTEST( basicStuff(Matrix()) ); CALL_SUBTEST( basicStuff(Matrix(10,10)) ); - + CALL_SUBTEST( basicStuffComplex(MatrixXcf(21, 17)) ); CALL_SUBTEST( basicStuffComplex(MatrixXcd(2, 3)) ); } From ab17f9272898390efff279c864a40f928752be43 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 10 Jul 2009 16:27:01 +0200 Subject: [PATCH 21/25] more sun studio fixes --- Eigen/src/Core/CwiseUnaryView.h | 2 +- Eigen/src/Core/Product.h | 6 ++++-- Eigen/src/Core/util/Constants.h | 7 ------- Eigen/src/Core/util/XprHelper.h | 11 +++++++++-- 4 files changed, 14 insertions(+), 12 deletions(-) diff --git a/Eigen/src/Core/CwiseUnaryView.h b/Eigen/src/Core/CwiseUnaryView.h index 8d8b51b4c..4785f01f4 100644 --- a/Eigen/src/Core/CwiseUnaryView.h +++ b/Eigen/src/Core/CwiseUnaryView.h @@ -45,7 +45,7 @@ struct ei_traits > ViewOp(typename ei_traits::Scalar) >::type Scalar; typedef typename MatrixType::Nested MatrixTypeNested; - typedef typename ei_unref::type _MatrixTypeNested; + typedef typename ei_cleantype::type _MatrixTypeNested; enum { Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), CoeffReadCost = ei_traits<_MatrixTypeNested>::CoeffReadCost + ei_functor_traits::Cost diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 2589d548e..c2317e425 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -83,13 +83,15 @@ struct ProductReturnType template struct ei_product_mode { enum{ + // workaround sun studio: + LhsIsVectorAtCompileTime = ei_traits::ColsAtCompileTime==1 || ei_traits::ColsAtCompileTime==1, value = ei_is_diagonal::ret || ei_is_diagonal::ret ? DiagonalProduct : ei_traits::MaxColsAtCompileTime == Dynamic && ( ei_traits::MaxRowsAtCompileTime == Dynamic || ei_traits::MaxColsAtCompileTime == Dynamic ) - && (!(ei_traits::IsVectorAtCompileTime && (ei_traits::Flags&RowMajorBit) && (!(ei_traits::Flags&DirectAccessBit)))) - && (!(ei_traits::IsVectorAtCompileTime && (!(ei_traits::Flags&RowMajorBit)) && (!(ei_traits::Flags&DirectAccessBit)))) + && (!(Rhs::IsVectorAtCompileTime && (ei_traits::Flags&RowMajorBit) && (!(ei_traits::Flags&DirectAccessBit)))) + && (!(LhsIsVectorAtCompileTime && (!(ei_traits::Flags&RowMajorBit)) && (!(ei_traits::Flags&DirectAccessBit)))) && (ei_is_same_type::Scalar, typename ei_traits::Scalar>::ret) ? CacheFriendlyProduct : NormalProduct }; diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 572cf7caf..585415662 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -198,13 +198,6 @@ const unsigned int SelfAdjoint = SelfAdjointBit; const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit; const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit; -template struct ei_is_diagonal -{ - enum { - ret = ( int(ei_traits::Flags) & DiagonalBits ) == DiagonalBits - }; -}; - enum { Aligned, Unaligned }; enum { ForceAligned, AsRequested }; enum { ConditionalJumpCost = 5 }; diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index b8d6aeb35..d8bbe11d5 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -217,7 +217,7 @@ template::ret > struct ei_special_scalar_op_base { - // dummy operator* so that the + // dummy operator* so that the // "using ei_special_scalar_op_base::operator*" compiles void operator*() const; }; @@ -249,7 +249,7 @@ template struct HNormalizedReturnType { typedef Block::ColsAtCompileTime==1 ? SizeMinusOne : 1, ei_traits::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne; - typedef CwiseUnaryOp::Scalar>, + typedef CwiseUnaryOp::Scalar>, NestByValue > Type; }; @@ -262,4 +262,11 @@ template struct ei_cast_return_type const XprType&,CastType>::ret type; }; +template struct ei_is_diagonal +{ + enum { + ret = ( int(ei_traits::Flags) & DiagonalBits ) == DiagonalBits + }; +}; + #endif // EIGEN_XPRHELPER_H From 1e7b1a8a85bf0ea54d4471c8b79de443f6bb0ccb Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 13 Jul 2009 14:55:03 +0200 Subject: [PATCH 22/25] add a SparseNestByValue expression and fix issue in sparse adjoint evaluation --- Eigen/Sparse | 1 + Eigen/src/Sparse/SparseMatrix.h | 5 +- Eigen/src/Sparse/SparseMatrixBase.h | 10 ++-- Eigen/src/Sparse/SparseNestByValue.h | 84 ++++++++++++++++++++++++++++ Eigen/src/Sparse/SparseUtil.h | 1 + test/main.h | 1 + test/sparse_basic.cpp | 2 + 7 files changed, 97 insertions(+), 7 deletions(-) create mode 100644 Eigen/src/Sparse/SparseNestByValue.h diff --git a/Eigen/Sparse b/Eigen/Sparse index 364fd50c9..a8888daa3 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -84,6 +84,7 @@ namespace Eigen { #include "src/Sparse/SparseUtil.h" #include "src/Sparse/SparseMatrixBase.h" +#include "src/Sparse/SparseNestByValue.h" #include "src/Sparse/CompressedStorage.h" #include "src/Sparse/AmbiVector.h" #include "src/Sparse/RandomSetter.h" diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index b751b5cb9..af3b5e5eb 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -443,9 +443,8 @@ class SparseMatrix // two passes algorithm: // 1 - compute the number of coeffs per dest inner vector // 2 - do the actual copy/eval - // Since each coeff of the rhs has to be evaluated twice, let's evauluate it if needed - //typedef typename ei_nested::type OtherCopy; - typedef typename ei_eval::type OtherCopy; + // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed + typedef typename ei_nested::type OtherCopy; typedef typename ei_cleantype::type _OtherCopy; OtherCopy otherCopy(other.derived()); diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 8d9406cfb..6cf4d5e96 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -99,8 +99,10 @@ template class SparseMatrixBase /** \internal the return type of MatrixBase::imag() */ typedef SparseCwiseUnaryOp, Derived> ImagReturnType; /** \internal the return type of MatrixBase::adjoint() */ - typedef SparseTranspose::type> /*>*/ - AdjointReturnType; + typedef typename ei_meta_if::IsComplex, + SparseCwiseUnaryOp, SparseNestByValue > >, + SparseTranspose + >::ret AdjointReturnType; #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is the "real scalar" type; if the \a Scalar type is already real numbers @@ -342,7 +344,7 @@ template class SparseMatrixBase SparseTranspose transpose() { return derived(); } const SparseTranspose transpose() const { return derived(); } // void transposeInPlace(); - const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; } + const AdjointReturnType adjoint() const { return transpose().nestByValue(); } // sub-vector SparseInnerVectorSet row(int i); @@ -520,7 +522,7 @@ template class SparseMatrixBase */ // inline int stride(void) const { return derived().stride(); } -// inline const NestByValue nestByValue() const; + inline const SparseNestByValue nestByValue() const; ConjugateReturnType conjugate() const; diff --git a/Eigen/src/Sparse/SparseNestByValue.h b/Eigen/src/Sparse/SparseNestByValue.h new file mode 100644 index 000000000..72fd056c4 --- /dev/null +++ b/Eigen/src/Sparse/SparseNestByValue.h @@ -0,0 +1,84 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2006-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_SPARSENESTBYVALUE_H +#define EIGEN_SPARSENESTBYVALUE_H + +/** \class SparseNestByValue + * + * \brief Expression which must be nested by value + * + * \param ExpressionType the type of the object of which we are requiring nesting-by-value + * + * This class is the return type of MatrixBase::nestByValue() + * and most of the time this is the only way it is used. + * + * \sa SparseMatrixBase::nestByValue(), class NestByValue + */ +template +struct ei_traits > : public ei_traits +{}; + +template class SparseNestByValue + : public SparseMatrixBase > +{ + public: + + class InnerIterator; + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseNestByValue) + + inline SparseNestByValue(const ExpressionType& matrix) : m_expression(matrix) {} + + EIGEN_STRONG_INLINE int rows() const { return m_expression.rows(); } + EIGEN_STRONG_INLINE int cols() const { return m_expression.cols(); } + + operator const ExpressionType&() const { return m_expression; } + + protected: + const ExpressionType m_expression; +}; + +/** \returns an expression of the temporary version of *this. + */ +template +inline const SparseNestByValue +SparseMatrixBase::nestByValue() const +{ + return SparseNestByValue(derived()); +} + +template +class SparseNestByValue::InnerIterator : public MatrixType::InnerIterator +{ + typedef typename MatrixType::InnerIterator Base; + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseNestByValue& expr, int outer) + : Base(expr.m_expression, outer) + {} +}; + +#endif // EIGEN_SPARSENESTBYVALUE_H diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index 52781aa46..a3b1d7ae9 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -106,6 +106,7 @@ template class DynamicSparseMatrix; template class SparseVector; template class MappedSparseMatrix; +template class SparseNestByValue; template class SparseTranspose; template class SparseInnerVectorSet; template class SparseCwise; diff --git a/test/main.h b/test/main.h index 53c8245c6..3947451cc 100644 --- a/test/main.h +++ b/test/main.h @@ -28,6 +28,7 @@ #include #include #include +#include #ifndef EIGEN_TEST_FUNC #error EIGEN_TEST_FUNC must be defined diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 9b7d7c4e6..666eea872 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -299,6 +299,8 @@ template void sparse_basic(const SparseMatrixType& re initSparse(density, refMat2, m2); VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); + + VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); } // test prune From bd506d837cf0ffa61352d72f4ca7a7c14bede9fa Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 13 Jul 2009 15:21:32 +0200 Subject: [PATCH 23/25] fix typo in previous commit --- Eigen/src/Sparse/SparseMatrixBase.h | 9 +-------- Eigen/src/Sparse/SparseNestByValue.h | 24 ++++++++++++------------ 2 files changed, 13 insertions(+), 20 deletions(-) diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 6cf4d5e96..83e79788c 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -446,14 +446,7 @@ template class SparseMatrixBase { for (typename Derived::InnerIterator i(derived(),j); i; ++i) { - if(IsRowMajor) - std::cerr << i.row() << "," << i.col() << " == " << j << "," << i.index() << "\n"; - else - std::cerr << i.row() << "," << i.col() << " == " << i.index() << "," << j << "\n"; -// if(IsRowMajor) - res.coeffRef(i.row(),i.col()) = i.value(); -// else -// res.coeffRef(i.index(),j) = i.value(); + res.coeffRef(i.row(),i.col()) = i.value(); } } return res; diff --git a/Eigen/src/Sparse/SparseNestByValue.h b/Eigen/src/Sparse/SparseNestByValue.h index 72fd056c4..b48277232 100644 --- a/Eigen/src/Sparse/SparseNestByValue.h +++ b/Eigen/src/Sparse/SparseNestByValue.h @@ -42,11 +42,11 @@ struct ei_traits > : public ei_traits class SparseNestByValue - : public SparseMatrixBase > + : public SparseMatrixBase > { public: - class InnerIterator; + typedef typename ExpressionType::InnerIterator InnerIterator; EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseNestByValue) @@ -70,15 +70,15 @@ SparseMatrixBase::nestByValue() const return SparseNestByValue(derived()); } -template -class SparseNestByValue::InnerIterator : public MatrixType::InnerIterator -{ - typedef typename MatrixType::InnerIterator Base; - public: - - EIGEN_STRONG_INLINE InnerIterator(const SparseNestByValue& expr, int outer) - : Base(expr.m_expression, outer) - {} -}; +// template +// class SparseNestByValue::InnerIterator : public MatrixType::InnerIterator +// { +// typedef typename MatrixType::InnerIterator Base; +// public: +// +// EIGEN_STRONG_INLINE InnerIterator(const SparseNestByValue& expr, int outer) +// : Base(expr.m_expression, outer) +// {} +// }; #endif // EIGEN_SPARSENESTBYVALUE_H From ddbaaebf9ee7bd1b6c3bb267ba5a1f3d6b63914a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 13 Jul 2009 15:27:01 +0200 Subject: [PATCH 24/25] one more fix of the previous commit (forgot to update ei_must_nest_by_value) --- Eigen/src/Sparse/SparseUtil.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index a3b1d7ae9..b5fc7c7b7 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -147,4 +147,6 @@ template class ei_eval typedef SparseMatrix<_Scalar, _Flags> type; }; +template struct ei_must_nest_by_value > { enum { ret = true }; }; + #endif // EIGEN_SPARSEUTIL_H From f5d2317b12f3d6eea68db5fb48743ca1801d10d3 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 13 Jul 2009 21:14:47 +0200 Subject: [PATCH 25/25] add a blueNorm() function implementing the Blues's stable norm algorithm. it is currently provided for experimentation purpose only. --- Eigen/src/Core/Dot.h | 138 +++++++++++++++++++++++++++++++++++- Eigen/src/Core/MatrixBase.h | 5 +- Eigen/src/Core/NumTraits.h | 8 ++- test/adjoint.cpp | 33 ++++++--- 4 files changed, 170 insertions(+), 14 deletions(-) diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 3861984da..4f185ea5b 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -295,7 +295,7 @@ inline typename NumTraits::Scalar>::Real MatrixBase< /** \returns the \em l2 norm of \c *this using a numerically more stable * algorithm. * - * \sa norm(), dot(), squaredNorm() + * \sa norm(), dot(), squaredNorm(), blueNorm() */ template inline typename NumTraits::Scalar>::Real @@ -304,6 +304,142 @@ MatrixBase::stableNorm() const return this->cwise().abs().redux(ei_scalar_hypot_op()); } +/** \internal Computes ibeta^iexp by binary expansion of iexp, + * exact if ibeta is the machine base */ +template inline T bexp(int ibeta, int iexp) +{ + T tbeta = T(ibeta); + T res = 1.0; + int n = iexp; + if (n<0) + { + n = - n; + tbeta = 1.0/tbeta; + } + for(;;) + { + if ((n % 2)==0) + res = res * tbeta; + n = n/2; + if (n==0) return res; + tbeta = tbeta*tbeta; + } + return res; +} + +/** \returns the \em l2 norm of \c *this using the Blue's algorithm. + * A Portable Fortran Program to Find the Euclidean Norm of a Vector, + * ACM TOMS, Vol 4, Issue 1, 1978. + * + * \sa norm(), dot(), squaredNorm(), stableNorm() + */ +template +inline typename NumTraits::Scalar>::Real +MatrixBase::blueNorm() const +{ + static int nmax; + static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr; + int n; + Scalar ax, abig, amed, asml; + + if(nmax <= 0) + { + int nbig, ibeta, it, iemin, iemax, iexp; + Scalar abig, eps; + // This program calculates the machine-dependent constants + // bl, b2, slm, s2m, relerr overfl, nmax + // from the "basic" machine-dependent numbers + // nbig, ibeta, it, iemin, iemax, rbig. + // The following define the basic machine-dependent constants. + // For portability, the PORT subprograms "ilmaeh" and "rlmach" + // are used. For any specific computer, each of the assignment + // statements can be replaced + nbig = std::numeric_limits::max(); // largest integer + ibeta = NumTraits::Base; // base for floating-point numbers + it = NumTraits::Mantissa; // number of base-beta digits in mantissa + iemin = std::numeric_limits::min_exponent; // minimum exponent + iemax = std::numeric_limits::max_exponent; // maximum exponent + rbig = std::numeric_limits::max(); // largest floating-point number + + // Check the basic machine-dependent constants. + if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) + || (it<=4 && ibeta <= 3 ) || it<2) + { + ei_assert(false && "the algorithm cannot be guaranteed on this computer"); + } + iexp = -((1-iemin)/2); + b1 = bexp(ibeta, iexp); // lower boundary of midrange + iexp = (iemax + 1 - it)/2; + b2 = bexp(ibeta,iexp); // upper boundary of midrange + + iexp = (2-iemin)/2; + s1m = bexp(ibeta,iexp); // scaling factor for lower range + iexp = - ((iemax+it)/2); + s2m = bexp(ibeta,iexp); // scaling factor for upper range + + overfl = rbig*s2m; // overfow boundary for abig + eps = bexp(ibeta, 1-it); + relerr = ei_sqrt(eps); // tolerance for neglecting asml + abig = 1.0/eps - 1.0; + if (Scalar(nbig)>abig) nmax = abig; // largest safe n + else nmax = nbig; + } + n = size(); + if(n==0) + return 0; + asml = Scalar(0); + amed = Scalar(0); + abig = Scalar(0); + for(int j=0; j b2) abig += ei_abs2(ax*s2m); + else if(ax < b2) asml += ei_abs2(ax*s1m); + else amed += ei_abs2(ax); + } + if(abig > Scalar(0)) + { + abig = ei_sqrt(abig); + if(abig > overfl) + { + ei_assert(false && "overflow"); + return rbig; + } + if(amed > Scalar(0)) + { + abig = abig/s2m; + amed = ei_sqrt(amed); + } + else + { + return abig/s2m; + } + + } + else if(asml > Scalar(0)) + { + if (amed > Scalar(0)) + { + abig = ei_sqrt(amed); + amed = ei_sqrt(asml) / s1m; + } + else + { + return ei_sqrt(asml)/s1m; + } + } + else + { + return ei_sqrt(amed); + } + asml = std::min(abig, amed); + abig = std::max(abig, amed); + if(asml <= abig*relerr) + return abig; + else + return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig)); +} + /** \returns an expression of the quotient of *this by its own norm. * * \only_for_vectors diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 2a49ae620..b8273ca22 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -381,8 +381,9 @@ template class MatrixBase template Scalar dot(const MatrixBase& other) const; RealScalar squaredNorm() const; - RealScalar norm() const; - RealScalar stableNorm() const; + RealScalar norm() const; + RealScalar stableNorm() const; + RealScalar blueNorm() const; const PlainMatrixType normalized() const; void normalize(); diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h index 24afe54c5..dec07a036 100644 --- a/Eigen/src/Core/NumTraits.h +++ b/Eigen/src/Core/NumTraits.h @@ -70,7 +70,9 @@ template<> struct NumTraits HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1 + MulCost = 1, + Base = 2, + Mantissa = 23 }; }; @@ -83,7 +85,9 @@ template<> struct NumTraits HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1 + MulCost = 1, + Base = 2, + Mantissa = 52 }; }; diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 965e4d256..14ee44a0f 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -76,6 +76,7 @@ template void adjoint(const MatrixType& m) { VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast(1)); VERIFY_IS_APPROX(v1.norm(), v1.stableNorm()); + VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); } // check compatibility of dot and adjoint @@ -113,15 +114,29 @@ template void adjoint(const MatrixType& m) void test_adjoint() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( adjoint(Matrix()) ); - CALL_SUBTEST( adjoint(Matrix3d()) ); - CALL_SUBTEST( adjoint(Matrix4f()) ); - CALL_SUBTEST( adjoint(MatrixXcf(4, 4)) ); - CALL_SUBTEST( adjoint(MatrixXi(8, 12)) ); - CALL_SUBTEST( adjoint(MatrixXf(21, 21)) ); - } +// for(int i = 0; i < g_repeat; i++) { +// CALL_SUBTEST( adjoint(Matrix()) ); +// CALL_SUBTEST( adjoint(Matrix3d()) ); +// CALL_SUBTEST( adjoint(Matrix4f()) ); +// CALL_SUBTEST( adjoint(MatrixXcf(4, 4)) ); +// CALL_SUBTEST( adjoint(MatrixXi(8, 12)) ); +// CALL_SUBTEST( adjoint(MatrixXf(21, 21)) ); +// } // test a large matrix only once - CALL_SUBTEST( adjoint(Matrix()) ); +// CALL_SUBTEST( adjoint(Matrix()) ); + for(int i = 0; i < g_repeat; i++) + { + std::cerr.precision(20); + int s = 1000000; + double y = 1.131242353467546478463457843445677435233e23 * ei_abs(ei_random()); + VectorXf v = VectorXf::Ones(s) * y; +// Vector4f x(v.segment(0,s/4).blueNorm(), v.segment(s/4+1,s/4).blueNorm(), +// v.segment((s/2)+1,s/4).blueNorm(), v.segment(3*s/4+1,s - 3*s/4-1).blueNorm()); +// std::cerr << v.norm() << " == " << v.stableNorm() << " == " << v.blueNorm() << " == " << x.norm() << "\n"; + std::cerr << v.norm() << "\n" << v.stableNorm() << "\n" << v.blueNorm() << "\n" << ei_sqrt(double(s)) * y << "\n\n\n"; + +// VectorXd d = VectorXd::Ones(s) * y;//v.cast(); +// std::cerr << d.norm() << "\n" << d.stableNorm() << "\n" << d.blueNorm() << "\n" << ei_sqrt(double(s)) * y << "\n\n\n"; + } }