From c7ae261ac08453e040caafb9383320c9ff95c48f Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 20 Aug 2009 01:29:38 -0400 Subject: [PATCH 01/16] adapt to API changes --- test/main.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/main.h b/test/main.h index 9ad014b64..aa563252f 100644 --- a/test/main.h +++ b/test/main.h @@ -227,10 +227,9 @@ inline bool test_ei_isMuchSmallerThan(const MatrixBase& m, return m.isMuchSmallerThan(s, test_precision::Scalar>()); } -template -void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::MatrixBase& m) +template +void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m) { - typedef Derived MatrixType; typedef typename ei_traits::Scalar Scalar; typedef Matrix VectorType; @@ -244,7 +243,7 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri HouseholderQR qra(a); HouseholderQR qrb(b); - m = (qra.matrixQ() * d * qrb.matrixQ()).lazy(); + m = qra.matrixQ() * d * qrb.matrixQ(); } } // end namespace Eigen From 72b002eab91cf058dfc667911db5e06cc1448a39 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 20 Aug 2009 12:19:15 -0400 Subject: [PATCH 02/16] work around internal compiler error with gcc 4.1 and 4.2, reported on the forum --- Eigen/src/Core/Product.h | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index e1e106b80..d70344deb 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -61,11 +61,22 @@ template struct ei_product_type enum { Rows = Lhs::RowsAtCompileTime, Cols = Rhs::ColsAtCompileTime, - Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime), - - value = ei_product_type_selector<(Rows >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Rows==1 ? 1 : Small)), - (Cols >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Cols==1 ? 1 : Small)), - (Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small))>::ret + Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime) + }; + + // the splitting into different lines of code here, introducing the _select enums and the typedef below, + // is to work around an internal compiler error with gcc 4.1 and 4.2. +private: + enum { + rows_select = Rows >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Rows==1 ? 1 : Small), + cols_select = Cols >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Cols==1 ? 1 : Small), + depth_select = Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small) + }; + typedef ei_product_type_selector product_type_selector; + +public: + enum { + value = product_type_selector::ret }; }; From 2f0b4e1abc6b105434fcb85df2461a37a810172b Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 21 Aug 2009 12:16:37 -0400 Subject: [PATCH 03/16] fix compilation with gcc 4.1. Indeed the path for recent gcc doesn't work with gcc 4.1, and looking at the implementation of vector in g++ 4.1, it was exactly our fallback case, so use that. --- Eigen/StdVector | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Eigen/StdVector b/Eigen/StdVector index a7f39ea37..f37de5ff6 100644 --- a/Eigen/StdVector +++ b/Eigen/StdVector @@ -136,10 +136,8 @@ class vector > { return vector_base::insert(position,x); } void insert(const_iterator position, size_type new_size, const value_type& x) { vector_base::insert(position, new_size, x); } -#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,1) +#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,2) // workaround GCC std::vector implementation - // Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&), - // no no need to workaround ! void resize(size_type new_size, const value_type& x) { if (new_size < vector_base::size()) @@ -147,9 +145,12 @@ class vector > else vector_base::insert(vector_base::end(), new_size - vector_base::size(), x); } -#elif defined(_GLIBCXX_VECTOR) +#elif defined(_GLIBCXX_VECTOR) && (!EIGEN_GNUC_AT_LEAST(4,1)) + // Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&), + // no no need to workaround ! using vector_base::resize; #else + // either GCC 4.1 or non-GCC // default implementation which should always work. void resize(size_type new_size, const value_type& x) { From ef582933c102fe514a161ccce950b77b403e7f2c Mon Sep 17 00:00:00 2001 From: "Marcus D. Hanwell" Date: Fri, 21 Aug 2009 14:04:17 -0400 Subject: [PATCH 04/16] Proper fix for linking to the Qt libraries (and others) My initial fix was incorrect, the libraries must be quoted when being passed to the add test macro, but must be unquoted when passed to the target_link_libraries function. --- cmake/EigenTesting.cmake | 2 +- test/CMakeLists.txt | 6 ------ 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake index ef36d2067..571774787 100644 --- a/cmake/EigenTesting.cmake +++ b/cmake/EigenTesting.cmake @@ -74,7 +74,7 @@ macro(ei_add_test testname) string(STRIP "${ARGV2}" ARGV2_stripped) string(LENGTH "${ARGV2_stripped}" ARGV2_stripped_length) if(${ARGV2_stripped_length} GREATER 0) - target_link_libraries(${targetname} "${ARGV2}") + target_link_libraries(${targetname} ${ARGV2}) endif(${ARGV2_stripped_length} GREATER 0) endif(${ARGC} GREATER 2) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ac4eb7633..dd8f07386 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -78,7 +78,6 @@ else(QT4_FOUND) ei_add_property(EIGEN_MISSING_BACKENDS "Qt4 support, ") endif(QT4_FOUND) - if(TEST_LIB) add_definitions("-DEIGEN_EXTERN_INSTANTIATIONS=1") endif(TEST_LIB) @@ -137,11 +136,6 @@ ei_add_test(regression) ei_add_test(stdvector) ei_add_test(resize) if(QT4_FOUND) - if(QT_QTCORE_LIBRARY_DEBUG) - set(QT_QTCORE_LIBRARY ${QT_QTCORE_LIBRARY_DEBUG}) - else(QT_QTCORE_LIBRARY_DEBUG) - set(QT_QTCORE_LIBRARY ${QT_QTCORE_LIBRARY_RELEASE}) - endif(QT_QTCORE_LIBRARY_DEBUG) ei_add_test(qtvector " " "${QT_QTCORE_LIBRARY}") endif(QT4_FOUND) ei_add_test(sparse_vector) From 7bedf5e9cb61a10e2be18a3f18ee63d68dbb4acd Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sat, 22 Aug 2009 01:13:21 -0400 Subject: [PATCH 05/16] add initial, rough, full-pivoting RRQR decomposition lots of room for improvement! and add Gael a (c) line in Householder.h --- Eigen/QR | 1 + Eigen/src/Householder/Householder.h | 1 + Eigen/src/QR/QR.h | 9 +- Eigen/src/QR/RRQR.h | 286 ++++++++++++++++++++++++++++ test/CMakeLists.txt | 1 + test/rrqr.cpp | 116 +++++++++++ 6 files changed, 410 insertions(+), 4 deletions(-) create mode 100644 Eigen/src/QR/RRQR.h create mode 100644 test/rrqr.cpp diff --git a/Eigen/QR b/Eigen/QR index 151c0b31b..c95f7522a 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -36,6 +36,7 @@ namespace Eigen { */ #include "src/QR/QR.h" +#include "src/QR/RRQR.h" #include "src/QR/Tridiagonalization.h" #include "src/QR/EigenSolver.h" #include "src/QR/SelfAdjointEigenSolver.h" diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 769ba3d43..8a274d240 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2009 Benoit Jacob +// Copyright (C) 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/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 90e6f8132..e5da6d691 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.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 @@ -39,7 +39,7 @@ * * Note that no pivoting is performed. This is \b not a rank-revealing decomposition. * - * \sa MatrixBase::qr() + * \sa MatrixBase::householderQr() */ template class HouseholderQR { @@ -54,6 +54,7 @@ template class HouseholderQR typedef Block MatrixRBlockType; typedef Matrix MatrixTypeR; typedef Matrix HCoeffsType; + typedef Matrix RowVectorType; /** * \brief Default Constructor. @@ -125,12 +126,12 @@ HouseholderQR& HouseholderQR::compute(const MatrixType& m_qr = matrix; m_hCoeffs.resize(size); - Matrix temp(cols); + RowVectorType temp(cols); for (int k = 0; k < size; ++k) { int remainingRows = rows - k; - int remainingCols = cols - k -1; + int remainingCols = cols - k - 1; RealScalar beta; m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); diff --git a/Eigen/src/QR/RRQR.h b/Eigen/src/QR/RRQR.h new file mode 100644 index 000000000..5e4f009dd --- /dev/null +++ b/Eigen/src/QR/RRQR.h @@ -0,0 +1,286 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2009 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_RRQR_H +#define EIGEN_RRQR_H + +/** \ingroup QR_Module + * \nonstableyet + * + * \class HouseholderRRQR + * + * \brief Householder rank-revealing QR decomposition of a matrix + * + * \param MatrixType the type of the matrix of which we are computing the QR decomposition + * + * This class performs a rank-revealing QR decomposition using Householder transformations. + * + * This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal + * numerical stability. + * + * \sa MatrixBase::householderRrqr() + */ +template class HouseholderRRQR +{ + public: + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) + }; + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix MatrixQType; + typedef Matrix HCoeffsType; + typedef Matrix IntRowVectorType; + typedef Matrix IntColVectorType; + typedef Matrix RowVectorType; + typedef Matrix ColVectorType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via HouseholderRRQR::compute(const MatrixType&). + */ + HouseholderRRQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + + HouseholderRRQR(const MatrixType& matrix) + : m_qr(matrix.rows(), matrix.cols()), + m_hCoeffs(std::min(matrix.rows(),matrix.cols())), + m_isInitialized(false) + { + compute(matrix); + } + + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which + * *this is the QR decomposition, if any exists. + * + * \param b the right-hand-side of the equation to solve. + * + * \param result a pointer to the vector/matrix in which to store the solution, if any exists. + * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). + * If no solution exists, *result is left with undefined coefficients. + * + * \note The case where b is a matrix is not yet implemented. Also, this + * code is space inefficient. + * + * Example: \include HouseholderRRQR_solve.cpp + * Output: \verbinclude HouseholderRRQR_solve.out + */ + template + 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 + */ + const MatrixType& matrixQR() const { return m_qr; } + + HouseholderRRQR& compute(const MatrixType& matrix); + + const IntRowVectorType& colsPermutation() const + { + ei_assert(m_isInitialized && "RRQR is not initialized."); + return m_cols_permutation; + } + + const IntColVectorType& rowsTranspositions() const + { + ei_assert(m_isInitialized && "RRQR is not initialized."); + return m_rows_transpositions; + } + + inline int rank() const + { + ei_assert(m_isInitialized && "RRQR is not initialized."); + return m_rank; + } + + protected: + MatrixType m_qr; + HCoeffsType m_hCoeffs; + IntColVectorType m_rows_transpositions; + IntRowVectorType m_cols_permutation; + bool m_isInitialized; + RealScalar m_precision; + int m_rank; + int m_det_pq; +}; + +#ifndef EIGEN_HIDE_HEAVY_CODE + +template +HouseholderRRQR& HouseholderRRQR::compute(const MatrixType& matrix) +{ + int rows = matrix.rows(); + int cols = matrix.cols(); + int size = std::min(rows,cols); + m_rank = size; + + m_qr = matrix; + m_hCoeffs.resize(size); + + RowVectorType temp(cols); + + // TODO: experiment to see the best formula + m_precision = epsilon() * size; + + m_rows_transpositions.resize(matrix.rows()); + IntRowVectorType cols_transpositions(matrix.cols()); + m_cols_permutation.resize(matrix.cols()); + int number_of_transpositions = 0; + + RealScalar biggest; + + for (int k = 0; k < size; ++k) + { + int row_of_biggest_in_corner, col_of_biggest_in_corner; + RealScalar biggest_in_corner; + + biggest_in_corner = m_qr.corner(Eigen::BottomRight, rows-k, cols-k) + .cwise().abs() + .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); + row_of_biggest_in_corner += k; + col_of_biggest_in_corner += k; + if(k==0) biggest = biggest_in_corner; + + // if the corner is negligible, then we have less than full rank, and we can finish early + if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) + { + m_rank = k; + for(int i = k; i < size; i++) + { + m_rows_transpositions.coeffRef(i) = i; + cols_transpositions.coeffRef(i) = i; + m_hCoeffs.coeffRef(i) = Scalar(0); + } + break; + } + + m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner; + cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; + if(k != row_of_biggest_in_corner) { + m_qr.row(k).end(cols-k).swap(m_qr.row(row_of_biggest_in_corner).end(cols-k)); + ++number_of_transpositions; + } + if(k != col_of_biggest_in_corner) { + m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner)); + ++number_of_transpositions; + } + + RealScalar beta; + m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); + m_qr.coeffRef(k,k) = beta; + + // apply H to remaining part of m_qr from the left + m_qr.corner(BottomRight, rows-k, cols-k-1) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + } + + for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; + for(int k = 0; k < size; ++k) + std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); + + m_det_pq = (number_of_transpositions%2) ? -1 : 1; + m_isInitialized = true; + + return *this; +} + +template +template +void HouseholderRRQR::solve( + const MatrixBase& b, + ResultType *result +) const +{ + ei_assert(m_isInitialized && "HouseholderRRQR is not initialized."); + const int rows = m_qr.rows(); + const int cols = b.cols(); + ei_assert(b.rows() == rows); + + typename OtherDerived::PlainMatrixType c(b); + + Matrix temp(cols); + for (int k = 0; k < m_rank; ++k) + { + int remainingSize = rows-k; + c.row(k).swap(c.row(m_rows_transpositions.coeff(k))); + c.corner(BottomRight, remainingSize, cols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); + } + + m_qr.corner(TopLeft, m_rank, m_rank) + .template triangularView() + .solveInPlace(c.corner(TopLeft, m_rank, c.cols())); + + result->resize(m_qr.cols(), b.cols()); + for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i); + for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero(); +} + +/** \returns the matrix Q */ +template +MatrixType HouseholderRRQR::matrixQ() const +{ + ei_assert(m_isInitialized && "HouseholderRRQR is not initialized."); + // compute the product H'_0 H'_1 ... H'_n-1, + // where H_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), ...] + int rows = m_qr.rows(); + int cols = m_qr.cols(); + int size = std::min(rows,cols); + MatrixType res = MatrixType::Identity(rows, rows); + Matrix temp(rows); + for (int k = size-1; k >= 0; k--) + { + res.block(k, k, rows-k, rows-k) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); + res.row(k).swap(res.row(m_rows_transpositions.coeff(k))); + } + return res; +} + +#endif // EIGEN_HIDE_HEAVY_CODE + +#if 0 +/** \return the Householder QR decomposition of \c *this. + * + * \sa class HouseholderRRQR + */ +template +const HouseholderRRQR::PlainMatrixType> +MatrixBase::householderQr() const +{ + return HouseholderRRQR(eval()); +} +#endif + + +#endif // EIGEN_QR_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index dd8f07386..b79e069b7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -121,6 +121,7 @@ ei_add_test(lu ${EI_OFLAG}) ei_add_test(determinant) ei_add_test(inverse ${EI_OFLAG}) ei_add_test(qr) +ei_add_test(rrqr) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") ei_add_test(svd) diff --git a/test/rrqr.cpp b/test/rrqr.cpp new file mode 100644 index 000000000..b6cc75d17 --- /dev/null +++ b/test/rrqr.cpp @@ -0,0 +1,116 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2009 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 . + +#include "main.h" +#include + +template void qr() +{ + /* this test covers the following files: QR.h */ + int rows = ei_random(20,200), cols = ei_random(20,200), cols2 = ei_random(20,200); + int rank = ei_random(1, std::min(rows, cols)-1); + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix SquareMatrixType; + typedef Matrix VectorType; + MatrixType m1; + createRandomMatrixOfRank(rank,rows,cols,m1); + HouseholderRRQR qr(m1); + VERIFY_IS_APPROX(rank, qr.rank()); + + MatrixType r = qr.matrixQR(); + // FIXME need better way to construct trapezoid + for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); + + MatrixType b = qr.matrixQ() * r; + + MatrixType c = MatrixType::Zero(rows,cols); + + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); + VERIFY_IS_APPROX(m1, c); + + MatrixType m2 = MatrixType::Random(cols,cols2); + MatrixType m3 = m1*m2; + m2 = MatrixType::Random(cols,cols2); + qr.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); +} + +template void qr_invertible() +{ + /* this test covers the following files: RRQR.h */ + typedef typename NumTraits::Real RealScalar; + int size = ei_random(10,50); + + MatrixType m1(size, size), m2(size, size), m3(size, size); + m1 = MatrixType::Random(size,size); + + if (ei_is_same_type::ret) + { + // let's build a matrix more stable to inverse + MatrixType a = MatrixType::Random(size,size*2); + m1 += a * a.adjoint(); + } + + HouseholderRRQR qr(m1); + m3 = MatrixType::Random(size,size); + qr.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); +} + +template void qr_verify_assert() +{ + MatrixType tmp; + + HouseholderRRQR qr; + VERIFY_RAISES_ASSERT(qr.matrixR()) + VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(qr.matrixQ()) +} + +void test_rrqr() +{ + for(int i = 0; i < 1; i++) { + // FIXME : very weird bug here +// CALL_SUBTEST( qr(Matrix2f()) ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + } + + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + } + + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); +} From 37dede6077250156be084e55e629c68535456d2a Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sat, 22 Aug 2009 10:40:39 -0400 Subject: [PATCH 06/16] fix typo --- test/stdvector.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/stdvector.cpp b/test/stdvector.cpp index 6779feae5..a976f0cf1 100644 --- a/test/stdvector.cpp +++ b/test/stdvector.cpp @@ -159,5 +159,5 @@ void test_stdvector() // some Quaternion CALL_SUBTEST(check_stdvector_quaternion(Quaternionf())); - CALL_SUBTEST(check_stdvector_quaternion(Quaternionf())); + CALL_SUBTEST(check_stdvector_quaternion(Quaterniond())); } From 90735b6a9cc4c91eb41ef1d2b6828f15a7f4fba3 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Sat, 22 Aug 2009 20:12:47 +0100 Subject: [PATCH 07/16] Rewrite tutorial section on solving linear systems --- doc/C01_QuickStartGuide.dox | 2 +- doc/C03_TutorialGeometry.dox | 2 +- doc/C05_TutorialLinearAlgebra.dox | 262 +++++++++++++----- doc/C07_TutorialSparse.dox | 2 +- doc/examples/Tutorial_PartialLU_solve.cpp | 18 ++ .../Tutorial_solve_matrix_inverse.cpp | 6 + doc/snippets/Tutorial_solve_multiple_rhs.cpp | 10 + .../Tutorial_solve_reuse_decomposition.cpp | 13 + doc/snippets/Tutorial_solve_singular.cpp | 9 + doc/snippets/Tutorial_solve_triangular.cpp | 8 + .../Tutorial_solve_triangular_inplace.cpp | 6 + 11 files changed, 267 insertions(+), 71 deletions(-) create mode 100644 doc/examples/Tutorial_PartialLU_solve.cpp create mode 100644 doc/snippets/Tutorial_solve_matrix_inverse.cpp create mode 100644 doc/snippets/Tutorial_solve_multiple_rhs.cpp create mode 100644 doc/snippets/Tutorial_solve_reuse_decomposition.cpp create mode 100644 doc/snippets/Tutorial_solve_singular.cpp create mode 100644 doc/snippets/Tutorial_solve_triangular.cpp create mode 100644 doc/snippets/Tutorial_solve_triangular_inplace.cpp diff --git a/doc/C01_QuickStartGuide.dox b/doc/C01_QuickStartGuide.dox index 3d98e14f5..d43dbd72d 100644 --- a/doc/C01_QuickStartGuide.dox +++ b/doc/C01_QuickStartGuide.dox @@ -1,6 +1,6 @@ namespace Eigen { -/** \page TutorialCore Tutorial 1/3 - Core features +/** \page TutorialCore Tutorial 1/4 - Core features \ingroup Tutorial
\ref index "Overview" diff --git a/doc/C03_TutorialGeometry.dox b/doc/C03_TutorialGeometry.dox index 2d3d742c9..e349c68ce 100644 --- a/doc/C03_TutorialGeometry.dox +++ b/doc/C03_TutorialGeometry.dox @@ -1,6 +1,6 @@ namespace Eigen { -/** \page TutorialGeometry Tutorial 2/3 - Geometry +/** \page TutorialGeometry Tutorial 2/4 - Geometry \ingroup Tutorial
\ref index "Overview" diff --git a/doc/C05_TutorialLinearAlgebra.dox b/doc/C05_TutorialLinearAlgebra.dox index 46333cb58..e70298b47 100644 --- a/doc/C05_TutorialLinearAlgebra.dox +++ b/doc/C05_TutorialLinearAlgebra.dox @@ -1,7 +1,6 @@ - namespace Eigen { -/** \page TutorialAdvancedLinearAlgebra Tutorial 3/3 - Advanced linear algebra +/** \page TutorialAdvancedLinearAlgebra Tutorial 3/4 - Advanced linear algebra \ingroup Tutorial
\ref index "Overview" @@ -11,6 +10,9 @@ namespace Eigen { | \ref TutorialSparse "Sparse matrix"
+This tutorial chapter explains how you can use Eigen to tackle various problems involving matrices: +solving systems of linear equations, finding eigenvalues and eigenvectors, and so on. + \b Table \b of \b contents - \ref TutorialAdvSolvers - \ref TutorialAdvLU @@ -18,53 +20,129 @@ namespace Eigen { - \ref TutorialAdvQR - \ref TutorialAdvEigenProblems + \section TutorialAdvSolvers Solving linear problems -This part of the tutorial focuses on solving linear problem of the form \f$ A \mathbf{x} = b \f$, -where both \f$ A \f$ and \f$ b \f$ are known, and \f$ x \f$ is the unknown. Moreover, \f$ A \f$ -assumed to be a squared matrix. Of course, the methods described here can be used whenever an expression -involve the product of an inverse matrix with a vector or another matrix: \f$ A^{-1} B \f$. -Eigen offers various algorithms to this problem, and its choice mainly depends on the nature of -the matrix \f$ A \f$, such as its shape, size and numerical properties. +This part of the tutorial focuses on solving systems of linear equations. Such statems can be +written in the form \f$ A \mathbf{x} = \mathbf{b} \f$, where both \f$ A \f$ and \f$ \mathbf{b} \f$ +are known, and \f$ \mathbf{x} \f$ is the unknown. Moreover, \f$ A \f$ is assumed to be a square +matrix. + +The equation \f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution if \f$ A \f$ is invertible. If +the matrix is not invertible, then the equation may have no or infinitely many solutions. All +solvers assume that \f$ A \f$ is invertible, unless noted otherwise. + +Eigen offers various algorithms to this problem. The choice of algorithm mainly depends on the +nature of the matrix \f$ A \f$, such as its shape, size and numerical properties. + - The \ref TutorialAdvSolvers_LU "LU decomposition" (with partial pivoting) is a general-purpose + algorithm which works for most problems. + - Use the \ref TutorialAdvSolvers_Cholesky "Cholesky decomposition" if the matrix \f$ A \f$ is + positive definite. + - Use a special \ref TutorialAdvSolvers_Triangular "triangular solver" if the matrix \f$ A \f$ is + upper or lower triangular. + - Use of the \ref TutorialAdvSolvers_Inverse "matrix inverse" is not recommended in general, but + may be appropriate in special cases, for instance if you want to solve several systems with the + same matrix \f$ A \f$ and that matrix is small. + - \ref TutorialAdvSolvers_Misc "Other solvers" (%LU decomposition with full pivoting, the singular + value decomposition) are provided for special cases, such as when \f$ A \f$ is not invertible. + +The methods described here can be used whenever an expression involve the product of an inverse +matrix with a vector or another matrix: \f$ A^{-1} \mathbf{v} \f$ or \f$ A^{-1} B \f$. + + +\subsection TutorialAdvSolvers_LU LU decomposition (with partial pivoting) + +This is a general-purpose algorithm which performs well in most cases (provided the matrix \f$ A \f$ +is invertible), so if you are unsure about which algorithm to pick, choose this. The method proceeds +in two steps. First, the %LU decomposition with partial pivoting is computed using the +MatrixBase::partialLu() function. This yields an object of the class PartialLU. Then, the +PartialLU::solve() method is called to compute a solution. + +As an example, suppose we want to solve the following system of linear equations: + +\f[ \begin{aligned} + x + 2y + 3z &= 3 \\ + 4x + 5y + 6z &= 3 \\ + 7x + 8y + 10z &= 4. +\end{aligned} \f] + +The following program solves this system: -\subsection TutorialAdvSolvers_Triangular Triangular solver -If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the diagonal -are all not zero), then the problem can be solved directly using MatrixBase::solveTriangular(), or better, -MatrixBase::solveTriangularInPlace(). Here is an example: -
-\include MatrixBase_marked.cpp - -output: -\include MatrixBase_marked.out +\include Tutorial_PartialLU_solve.cpp + +output: \include Tutorial_PartialLU_solve.out
-See MatrixBase::solveTriangular() for more details. +There are many situations in which we want to solve the same system of equations with different +right-hand sides. One possibility is to put the right-hand sides as columns in a matrix \f$ B \f$ +and then solve the equation \f$ A X = B \f$. For instance, suppose that we want to solve the same +system as before, but now we also need the solution of the same equations with 1 on the right-hand +side. The following code computes the required solutions: + +
+\include Tutorial_solve_multiple_rhs.cpp + +output: \include Tutorial_solve_multiple_rhs.out +
+ +However, this is not always possible. Often, you only know the right-hand side of the second +problem, and whether you want to solve it at all, after you solved the first problem. In such a +case, it's best to save the %LU decomposition and reuse it to solve the second problem. This is +worth the effort because computing the %LU decomposition is much more expensive than using it to +solve the equation. Here is some code to illustrate the procedure. It uses the constructor +PartialLU::PartialLU(const MatrixType&) to compute the %LU decomposition. + +
+\include Tutorial_solve_reuse_decomposition.cpp + +output: \include Tutorial_solve_reuse_decomposition.out +
+ +\b Warning: All this code presumes that the matrix \f$ A \f$ is invertible, so that the system +\f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution. If the matrix \f$ A \f$ is not invertible, +then the system \f$ A \mathbf{x} = \mathbf{b} \f$ has either zero or infinitely many solutions. In +both cases, PartialLU::solve() will give nonsense results. For example, suppose that we want to +solve the same system as above, but with the 10 in the last equation replaced by 9. Then the system +of equations is inconsistent: adding the first and the third equation gives \f$ 8x + 10y + 12z = 7 \f$, +which implies \f$ 4x + 5y + 6z = 3\frac12 \f$, in contradiction with the seocond equation. If we try +to solve this inconsistent system with Eigen, we find: + +
+\include Tutorial_solve_singular.cpp + +output: \include Tutorial_solve_singular.out +
+ +The %LU decomposition with \b full pivoting (class LU) and the singular value decomposition (class +SVD) may be helpful in this case, as explained in the section \ref TutorialAdvSolvers_Misc below. + +\sa LU_Module, MatrixBase::partialLu(), PartialLU::solve(), class PartialLU. -\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices) -If the matrix \f$ A \f$ is small (\f$ \leq 4 \f$) and invertible, then a good approach is to directly compute -the inverse of the matrix \f$ A \f$, and then obtain the solution \f$ x \f$ by \f$ \mathbf{x} = A^{-1} b \f$. With Eigen, -this can be implemented like this: +\subsection TutorialAdvSolvers_Cholesky Cholesky decomposition -\code -#include -Matrix4f A = Matrix4f::Random(); -Vector4f b = Vector4f::Random(); -Vector4f x = A.inverse() * b; -\endcode +If the matrix \f$ A \f$ is \b symmetric \b positive \b definite, then the best method is to use a +Cholesky decomposition: it is both faster and more accurate than the %LU decomposition. Such +positive definite matrices often arise when solving overdetermined problems. These are linear +systems \f$ A \mathbf{x} = \mathbf{b} \f$ in which the matrix \f$ A \f$ has more rows than columns, +so that there are more equations than unknowns. Typically, there is no vector \f$ \mathbf{x} \f$ +which satisfies all the equation. Instead, we look for the least-square solution, that is, the +vector \f$ \mathbf{x} \f$ for which \f$ \| A \mathbf{x} - \mathbf{b} \|_2 \f$ is minimal. You can +find this vector by solving the equation \f$ A^T \! A \mathbf{x} = A^T \mathbf{b} \f$. If the matrix +\f$ A \f$ has full rank, then \f$ A^T \! A \f$ is positive definite and thus you can use the +Cholesky decomposition to solve this system and find the least-square solution to the original +system \f$ A \mathbf{x} = \mathbf{b} \f$. -Note that the function inverse() is defined in the LU module. -See MatrixBase::inverse() for more details. +Eigen offers two different Cholesky decompositions: the LLT class provides a \f$ LL^T \f$ +decomposition where L is a lower triangular matrix, and the LDLT class provides a \f$ LDL^T \f$ +decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. The latter +includes pivoting and avoids square roots; this makes the %LDLT decomposition slightly more stable +than the %LLT decomposition. The LDLT class is able to handle both positive- and negative-definite +matrices, but not indefinite matrices. +The API is the same as when using the %LU decomposition. -\subsection TutorialAdvSolvers_Symmetric Cholesky (for positive definite matrices) -If the matrix \f$ A \f$ is \b positive \b definite, then -the best method is to use a Cholesky decomposition. -Such positive definite matrices often arise when solving overdetermined problems in a least square sense (see below). -Eigen offers two different Cholesky decompositions: a \f$ LL^T \f$ decomposition where L is a lower triangular matrix, -and a \f$ LDL^T \f$ decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. -The latter avoids square roots and is therefore slightly more stable than the former one. \code #include MatrixXf D = MatrixXf::Random(8,4); @@ -74,15 +152,19 @@ VectorXf x; A.llt().solve(b,&x); // using a LLT factorization A.ldlt().solve(b,&x); // using a LDLT factorization \endcode -when the value of the right hand side \f$ b \f$ is not needed anymore, then it is faster to use -the \em in \em place API, e.g.: + +The LLT and LDLT classes also provide an \em in \em place API for the case where the value of the +right hand-side \f$ b \f$ is not needed anymore. + \code A.llt().solveInPlace(b); \endcode -In that case the value of \f$ b \f$ is replaced by the result \f$ x \f$. -If the linear problem has to solved for various right hand side \f$ b_i \f$, but same matrix \f$ A \f$, -then it is highly recommended to perform the decomposition of \$ A \$ only once, e.g.: +This code replaces the vector \f$ b \f$ by the result \f$ x \f$. + +As before, you can reuse the factorization if you have to solve the same linear problem with +different right-hand sides, e.g.: + \code // ... LLT lltOfA(A); @@ -91,40 +173,69 @@ lltOfA.solveInPlace(b1); // ... \endcode -\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT. +\sa Cholesky_Module, MatrixBase::llt(), MatrixBase::ldlt(), LLT::solve(), LLT::solveInPlace(), +LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT. -\subsection TutorialAdvSolvers_LU LU decomposition (for most cases) -If the matrix \f$ A \f$ does not fit in any of the previous categories, or if you are unsure about the numerical -stability of your problem, then you can use the LU solver based on a decomposition of the same name : see the section \ref TutorialAdvLU below. Actually, Eigen's LU module does not implement a standard LU decomposition, but rather a so-called LU decomposition -with full pivoting and rank update which has much better numerical stability. -The API of the LU solver is the same than the Cholesky one, except that there is no \em in \em place variant: -\code -#include -MatrixXf A = MatrixXf::Random(20,20); -VectorXf b = VectorXf::Random(20); -VectorXf x; -A.lu().solve(b, &x); -\endcode +\subsection TutorialAdvSolvers_Triangular Triangular solver -Again, the LU decomposition can be stored to be reused or to perform other kernel operations: -\code -// ... -LU luOfA(A); -luOfA.solve(b, &x); -// ... -\endcode +If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the +diagonal are all not zero), then the problem can be solved directly using the TriangularView +class. This is much faster than using an %LU or Cholesky decomposition (in fact, the triangular +solver is used when you solve a system using the %LU or Cholesky decomposition). Here is an example: -See the section \ref TutorialAdvLU below. +
+\include Tutorial_solve_triangular.cpp + +output: \include Tutorial_solve_triangular.out +
-\sa class LU, LU::solve(), LU_Module +The MatrixBase::triangularView() function constructs an object of the class TriangularView, and +TriangularView::solve() then solves the system. There is also an \e in \e place variant: + +
+\include Tutorial_solve_triangular_inplace.cpp + +output: \include Tutorial_solve_triangular_inplace.out +
+ +\sa MatrixBase::triangularView(), TriangularView::solve(), TriangularView::solveInPlace(), +TriangularView class. -\subsection TutorialAdvSolvers_SVD SVD solver (for singular matrices and special cases) -Finally, Eigen also offer a solver based on a singular value decomposition (SVD). Again, the API is the -same than with Cholesky or LU: +\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices) + +The solution of the system \f$ A \mathbf{x} = \mathbf{b} \f$ is given by \f$ \mathbf{x} = A^{-1} +\mathbf{b} \f$. This suggests the following approach for solving the system: compute the matrix +inverse and multiply that with the right-hand side. This is often not a good approach: using the %LU +decomposition with partial pivoting yields a more accurate algorithm that is usually just as fast or +even faster. However, using the matrix inverse can be faster if the matrix \f$ A \f$ is small +(≤4) and fixed size, though numerical stability problems may still remain. Here is an example of +how you would write this in Eigen: + +
+\include Tutorial_solve_matrix_inverse.cpp + +output: \include Tutorial_solve_matrix_inverse.out +
+ +Note that the function inverse() is defined in the \ref LU_Module. + +\sa MatrixBase::inverse(). + + +\subsection TutorialAdvSolvers_Misc Other solvers (for singular matrices and special cases) + +Finally, Eigen also offer solvers based on a singular value decomposition (%SVD) or the %LU +decomposition with full pivoting. These have the same API as the solvers based on the %LU +decomposition with partial pivoting (PartialLU). + +The solver based on the %SVD uses the class SVD. It can handle singular matrices. Here is an example +of its use: + \code #include +// ... MatrixXf A = MatrixXf::Random(20,20); VectorXf b = VectorXf::Random(20); VectorXf x; @@ -133,8 +244,23 @@ SVD svdOfA(A); svdOfA.solve(b, &x); \endcode -\sa class SVD, SVD::solve(), SVD_Module +%LU decomposition with full pivoting has better numerical stability than %LU decomposition with +partial pivoting. It is defined in the class LU. The solver can also handle singular matrices. +\code +#include +// ... +MatrixXf A = MatrixXf::Random(20,20); +VectorXf b = VectorXf::Random(20); +VectorXf x; +A.lu().solve(b, &x); +LU luOfA(A); +luOfA.solve(b, &x); +\endcode + +See the section \ref TutorialAdvLU below. + +\sa class SVD, SVD::solve(), SVD_Module, class LU, LU::solve(), LU_Module. diff --git a/doc/C07_TutorialSparse.dox b/doc/C07_TutorialSparse.dox index a8bfe006e..3a7182883 100644 --- a/doc/C07_TutorialSparse.dox +++ b/doc/C07_TutorialSparse.dox @@ -1,6 +1,6 @@ namespace Eigen { -/** \page TutorialSparse Tutorial - Getting started with the sparse module +/** \page TutorialSparse Tutorial 4/4 - Getting started with the sparse module \ingroup Tutorial
\ref index "Overview" diff --git a/doc/examples/Tutorial_PartialLU_solve.cpp b/doc/examples/Tutorial_PartialLU_solve.cpp new file mode 100644 index 000000000..80c393f9a --- /dev/null +++ b/doc/examples/Tutorial_PartialLU_solve.cpp @@ -0,0 +1,18 @@ +#include +#include + +using namespace std; +using namespace Eigen; + +int main(int, char *[]) +{ + Matrix3f A; + Vector3f b; + A << 1,2,3, 4,5,6, 7,8,10; + b << 3, 3, 4; + cout << "Here is the matrix A:" << endl << A << endl; + cout << "Here is the vector b:" << endl << b << endl; + Vector3f x; + A.partialLu().solve(b, &x); + cout << "The solution is:" << endl << x << endl; +} diff --git a/doc/snippets/Tutorial_solve_matrix_inverse.cpp b/doc/snippets/Tutorial_solve_matrix_inverse.cpp new file mode 100644 index 000000000..fff324446 --- /dev/null +++ b/doc/snippets/Tutorial_solve_matrix_inverse.cpp @@ -0,0 +1,6 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 4,5,6, 7,8,10; +b << 3, 3, 4; +Vector3f x = A.inverse() * b; +cout << "The solution is:" << endl << x << endl; diff --git a/doc/snippets/Tutorial_solve_multiple_rhs.cpp b/doc/snippets/Tutorial_solve_multiple_rhs.cpp new file mode 100644 index 000000000..fbb15165a --- /dev/null +++ b/doc/snippets/Tutorial_solve_multiple_rhs.cpp @@ -0,0 +1,10 @@ +Matrix3f A(3,3); +A << 1,2,3, 4,5,6, 7,8,10; +Matrix B; +B << 3,1, 3,1, 4,1; +Matrix X; +A.partialLu().solve(B, &X); +cout << "The solution with right-hand side (3,3,4) is:" << endl; +cout << X.col(0) << endl; +cout << "The solution with right-hand side (1,1,1) is:" << endl; +cout << X.col(1) << endl; diff --git a/doc/snippets/Tutorial_solve_reuse_decomposition.cpp b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp new file mode 100644 index 000000000..b4112adc4 --- /dev/null +++ b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp @@ -0,0 +1,13 @@ +Matrix3f A(3,3); +A << 1,2,3, 4,5,6, 7,8,10; +PartialLU luOfA(A); // compute LU decomposition of A +Vector3f b; +b << 3,3,4; +Vector3f x; +luOfA.solve(b, &x); +cout << "The solution with right-hand side (3,3,4) is:" << endl; +cout << x << endl; +b << 1,1,1; +luOfA.solve(b, &x); +cout << "The solution with right-hand side (1,1,1) is:" << endl; +cout << x << endl; diff --git a/doc/snippets/Tutorial_solve_singular.cpp b/doc/snippets/Tutorial_solve_singular.cpp new file mode 100644 index 000000000..da94ad445 --- /dev/null +++ b/doc/snippets/Tutorial_solve_singular.cpp @@ -0,0 +1,9 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 4,5,6, 7,8,9; +b << 3, 3, 4; +cout << "Here is the matrix A:" << endl << A << endl; +cout << "Here is the vector b:" << endl << b << endl; +Vector3f x; +A.partialLu().solve(b, &x); +cout << "The solution is:" << endl << x << endl; diff --git a/doc/snippets/Tutorial_solve_triangular.cpp b/doc/snippets/Tutorial_solve_triangular.cpp new file mode 100644 index 000000000..ff876f687 --- /dev/null +++ b/doc/snippets/Tutorial_solve_triangular.cpp @@ -0,0 +1,8 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 0,5,6, 0,0,10; +b << 3, 3, 4; +cout << "Here is the matrix A:" << endl << A << endl; +cout << "Here is the vector b:" << endl << b << endl; +Vector3f x = A.triangularView().solve(b); +cout << "The solution is:" << endl << x << endl; diff --git a/doc/snippets/Tutorial_solve_triangular_inplace.cpp b/doc/snippets/Tutorial_solve_triangular_inplace.cpp new file mode 100644 index 000000000..64ae4eea1 --- /dev/null +++ b/doc/snippets/Tutorial_solve_triangular_inplace.cpp @@ -0,0 +1,6 @@ +Matrix3f A; +Vector3f b; +A << 1,2,3, 0,5,6, 0,0,10; +b << 3, 3, 4; +A.triangularView().solveInPlace(b); +cout << "The solution is:" << endl << b << endl; From a848ed02ad51f3223e11d7c39b59b9e2c16b634f Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 23 Aug 2009 17:33:31 -0400 Subject: [PATCH 08/16] let createRandomMatrixOfRank support fixed-size! --- test/main.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/main.h b/test/main.h index aa563252f..e3866c68b 100644 --- a/test/main.h +++ b/test/main.h @@ -231,7 +231,7 @@ template void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m) { typedef typename ei_traits::Scalar Scalar; - typedef Matrix VectorType; + typedef Matrix VectorType; MatrixType a = MatrixType::Random(rows,rows); MatrixType d = MatrixType::Identity(rows,cols); From e86dbd5255ade8c875ab2712667224e9745c62d9 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 23 Aug 2009 23:49:44 +0200 Subject: [PATCH 09/16] fix apply Jacobi for complexes and add documentation for some *Jacobi* functions --- Eigen/src/Jacobi/Jacobi.h | 36 +++++++++++++++++++++++++++++------- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index b5940c74b..28b6cc7ad 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -26,9 +26,20 @@ #ifndef EIGEN_JACOBI_H #define EIGEN_JACOBI_H +/** Applies the counter clock wise 2D rotation of angle \c theta given by its + * cosine \a c and sine \a s to the set of 2D vectors of cordinates \a x and \a y: + * \f$ x = c x - s' y \f$ + * \f$ y = s x + c y \f$ + * + * \sa MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() + */ template void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s); +/** Applies a rotation in the plane defined by \a c, \a s to the rows \a p and \a q of \c *this. + * More precisely, it computes B = J' * B, with J = [c s ; -s' c] and B = [ *this.row(p) ; *this.row(q) ] + * \sa MatrixBase::applyJacobiOnTheRight(), ei_apply_rotation_in_the_plane() + */ template inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) { @@ -37,6 +48,10 @@ inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Sc ei_apply_rotation_in_the_plane(x, y, c, s); } +/** Applies a rotation in the plane defined by \a c, \a s to the columns \a p and \a q of \c *this. + * More precisely, it computes B = B * J, with J = [c s ; -s' c] and B = [ *this.col(p) ; *this.col(q) ] + * \sa MatrixBase::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane() + */ template inline void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) { @@ -45,6 +60,12 @@ inline void MatrixBase::applyJacobiOnTheRight(int p, int q, Scalar c, S ei_apply_rotation_in_the_plane(x, y, c, s); } +/** Computes the cosine-sine pair (\a c, \a s) such that its associated + * rotation \f$ J = ( \begin{array}{cc} c & s \\ -s' c \end{array} )\f$ + * applied to both the right and left of the 2x2 matrix + * \f$ B = ( \begin{array}{cc} x & y \\ _ & z \end{array} )\f$ yields + * a diagonal matrix A: \f$ A = J' B J \f$ + */ template bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) { @@ -128,12 +149,13 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& const Packet pc = ei_pset1(c); const Packet ps = ei_pset1(s); + ei_conj_helper cj; for(int i=0; i Date: Mon, 24 Aug 2009 00:01:02 +0200 Subject: [PATCH 10/16] hm, forgot to conjugate the arguments in applyJacobiOnTheLeft --- Eigen/src/Jacobi/Jacobi.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 28b6cc7ad..a6670a49f 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -45,7 +45,7 @@ inline void MatrixBase::applyJacobiOnTheLeft(int p, int q, Scalar c, Sc { RowXpr x(row(p)); RowXpr y(row(q)); - ei_apply_rotation_in_the_plane(x, y, c, s); + ei_apply_rotation_in_the_plane(x, y, ei_conj(c), ei_conj(s)); } /** Applies a rotation in the plane defined by \a c, \a s to the columns \a p and \a q of \c *this. From 97bc1af1f107463e23de4b5207202ded478cddaa Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 23 Aug 2009 18:04:33 -0400 Subject: [PATCH 11/16] add ColPivotingHouseholderQR rename RRQR to fullPivotingHouseholderQR --- Eigen/QR | 3 +- Eigen/src/Core/MatrixBase.h | 2 + Eigen/src/Core/util/ForwardDeclarations.h | 2 + Eigen/src/QR/ColPivotingHouseholderQR.h | 266 ++++++++++++++++++ .../{RRQR.h => FullPivotingHouseholderQR.h} | 53 ++-- test/CMakeLists.txt | 3 +- test/{rrqr.cpp => qr_colpivoting.cpp} | 8 +- test/qr_fullpivoting.cpp | 115 ++++++++ 8 files changed, 417 insertions(+), 35 deletions(-) create mode 100644 Eigen/src/QR/ColPivotingHouseholderQR.h rename Eigen/src/QR/{RRQR.h => FullPivotingHouseholderQR.h} (84%) rename test/{rrqr.cpp => qr_colpivoting.cpp} (95%) create mode 100644 test/qr_fullpivoting.cpp diff --git a/Eigen/QR b/Eigen/QR index c95f7522a..a0575040c 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -36,7 +36,8 @@ namespace Eigen { */ #include "src/QR/QR.h" -#include "src/QR/RRQR.h" +#include "src/QR/FullPivotingHouseholderQR.h" +#include "src/QR/ColPivotingHouseholderQR.h" #include "src/QR/Tridiagonalization.h" #include "src/QR/EigenSolver.h" #include "src/QR/SelfAdjointEigenSolver.h" diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 1f4c6bf7a..7cb8853e6 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -745,6 +745,8 @@ template class MatrixBase /////////// QR module /////////// const HouseholderQR householderQr() const; + const ColPivotingHouseholderQR colPivotingHouseholderQr() const; + const FullPivotingHouseholderQR fullPivotingHouseholderQr() const; EigenvaluesReturnType eigenvalues() const; RealScalar operatorNorm() const; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 6f3e6a788..414aaceca 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -117,6 +117,8 @@ template class Reverse; template class LU; template class PartialLU; template class HouseholderQR; +template class ColPivotingHouseholderQR; +template class FullPivotingHouseholderQR; template class SVD; template class LLT; template class LDLT; diff --git a/Eigen/src/QR/ColPivotingHouseholderQR.h b/Eigen/src/QR/ColPivotingHouseholderQR.h new file mode 100644 index 000000000..ed4b84f63 --- /dev/null +++ b/Eigen/src/QR/ColPivotingHouseholderQR.h @@ -0,0 +1,266 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2009 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_COLPIVOTINGHOUSEHOLDERQR_H +#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H + +/** \ingroup QR_Module + * \nonstableyet + * + * \class ColPivotingHouseholderQR + * + * \brief Householder rank-revealing QR decomposition of a matrix + * + * \param MatrixType the type of the matrix of which we are computing the QR decomposition + * + * This class performs a rank-revealing QR decomposition using Householder transformations. + * + * This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal + * numerical stability. + * + * \sa MatrixBase::colPivotingHouseholderQr() + */ +template class ColPivotingHouseholderQR +{ + public: + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) + }; + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix MatrixQType; + typedef Matrix HCoeffsType; + typedef Matrix IntRowVectorType; + typedef Matrix IntColVectorType; + typedef Matrix RowVectorType; + typedef Matrix ColVectorType; + typedef Matrix RealRowVectorType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ColPivotingHouseholderQR::compute(const MatrixType&). + */ + ColPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + + ColPivotingHouseholderQR(const MatrixType& matrix) + : m_qr(matrix.rows(), matrix.cols()), + m_hCoeffs(std::min(matrix.rows(),matrix.cols())), + m_isInitialized(false) + { + compute(matrix); + } + + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which + * *this is the QR decomposition, if any exists. + * + * \param b the right-hand-side of the equation to solve. + * + * \param result a pointer to the vector/matrix in which to store the solution, if any exists. + * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). + * If no solution exists, *result is left with undefined coefficients. + * + * \note The case where b is a matrix is not yet implemented. Also, this + * code is space inefficient. + * + * Example: \include ColPivotingHouseholderQR_solve.cpp + * Output: \verbinclude ColPivotingHouseholderQR_solve.out + */ + template + 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 + */ + const MatrixType& matrixQR() const { return m_qr; } + + ColPivotingHouseholderQR& compute(const MatrixType& matrix); + + const IntRowVectorType& colsPermutation() const + { + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + return m_cols_permutation; + } + + inline int rank() const + { + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + return m_rank; + } + + protected: + MatrixType m_qr; + HCoeffsType m_hCoeffs; + IntRowVectorType m_cols_permutation; + bool m_isInitialized; + RealScalar m_precision; + int m_rank; + int m_det_pq; +}; + +#ifndef EIGEN_HIDE_HEAVY_CODE + +template +ColPivotingHouseholderQR& ColPivotingHouseholderQR::compute(const MatrixType& matrix) +{ + int rows = matrix.rows(); + int cols = matrix.cols(); + int size = std::min(rows,cols); + m_rank = size; + + m_qr = matrix; + m_hCoeffs.resize(size); + + RowVectorType temp(cols); + + m_precision = epsilon() * size; + + IntRowVectorType cols_transpositions(matrix.cols()); + m_cols_permutation.resize(matrix.cols()); + int number_of_transpositions = 0; + + RealRowVectorType colSqNorms(cols); + for(int k = 0; k < cols; ++k) + colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); + RealScalar biggestColSqNorm = colSqNorms.maxCoeff(); + + for (int k = 0; k < size; ++k) + { + int biggest_col_in_corner; + RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner); + biggest_col_in_corner += k; + + // if the corner is negligible, then we have less than full rank, and we can finish early + if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision)) + { + m_rank = k; + for(int i = k; i < size; i++) + { + cols_transpositions.coeffRef(i) = i; + m_hCoeffs.coeffRef(i) = Scalar(0); + } + break; + } + + cols_transpositions.coeffRef(k) = biggest_col_in_corner; + if(k != biggest_col_in_corner) { + m_qr.col(k).swap(m_qr.col(biggest_col_in_corner)); + ++number_of_transpositions; + } + + RealScalar beta; + m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); + m_qr.coeffRef(k,k) = beta; + + m_qr.corner(BottomRight, rows-k, cols-k-1) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + + colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2(); + } + + for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; + for(int k = 0; k < size; ++k) + std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); + + m_det_pq = (number_of_transpositions%2) ? -1 : 1; + m_isInitialized = true; + + return *this; +} + +template +template +void ColPivotingHouseholderQR::solve( + const MatrixBase& b, + ResultType *result +) const +{ + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + const int rows = m_qr.rows(); + const int cols = b.cols(); + ei_assert(b.rows() == rows); + + typename OtherDerived::PlainMatrixType c(b); + + Matrix temp(cols); + for (int k = 0; k < m_rank; ++k) + { + int remainingSize = rows-k; + c.corner(BottomRight, remainingSize, cols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); + } + + m_qr.corner(TopLeft, m_rank, m_rank) + .template triangularView() + .solveInPlace(c.corner(TopLeft, m_rank, c.cols())); + + result->resize(m_qr.cols(), b.cols()); + for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i); + for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero(); +} + +/** \returns the matrix Q */ +template +MatrixType ColPivotingHouseholderQR::matrixQ() const +{ + ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); + // compute the product H'_0 H'_1 ... H'_n-1, + // where H_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), ...] + int rows = m_qr.rows(); + int cols = m_qr.cols(); + int size = std::min(rows,cols); + MatrixType res = MatrixType::Identity(rows, rows); + Matrix temp(rows); + for (int k = size-1; k >= 0; k--) + { + res.block(k, k, rows-k, rows-k) + .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); + } + return res; +} + +#endif // EIGEN_HIDE_HEAVY_CODE + +/** \return the column-pivoting Householder QR decomposition of \c *this. + * + * \sa class ColPivotingHouseholderQR + */ +template +const ColPivotingHouseholderQR::PlainMatrixType> +MatrixBase::colPivotingHouseholderQr() const +{ + return ColPivotingHouseholderQR(eval()); +} + + +#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H diff --git a/Eigen/src/QR/RRQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h similarity index 84% rename from Eigen/src/QR/RRQR.h rename to Eigen/src/QR/FullPivotingHouseholderQR.h index 5e4f009dd..f7b0f1cc1 100644 --- a/Eigen/src/QR/RRQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -23,13 +23,13 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#ifndef EIGEN_RRQR_H -#define EIGEN_RRQR_H +#ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H +#define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H /** \ingroup QR_Module * \nonstableyet * - * \class HouseholderRRQR + * \class FullPivotingHouseholderQR * * \brief Householder rank-revealing QR decomposition of a matrix * @@ -42,7 +42,7 @@ * * \sa MatrixBase::householderRrqr() */ -template class HouseholderRRQR +template class FullPivotingHouseholderQR { public: @@ -66,11 +66,11 @@ template class HouseholderRRQR * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to - * perform decompositions via HouseholderRRQR::compute(const MatrixType&). + * perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&). */ - HouseholderRRQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + FullPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} - HouseholderRRQR(const MatrixType& matrix) + FullPivotingHouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs(std::min(matrix.rows(),matrix.cols())), m_isInitialized(false) @@ -90,8 +90,8 @@ template class HouseholderRRQR * \note The case where b is a matrix is not yet implemented. Also, this * code is space inefficient. * - * Example: \include HouseholderRRQR_solve.cpp - * Output: \verbinclude HouseholderRRQR_solve.out + * Example: \include FullPivotingHouseholderQR_solve.cpp + * Output: \verbinclude FullPivotingHouseholderQR_solve.out */ template void solve(const MatrixBase& b, ResultType *result) const; @@ -102,23 +102,23 @@ template class HouseholderRRQR */ const MatrixType& matrixQR() const { return m_qr; } - HouseholderRRQR& compute(const MatrixType& matrix); + FullPivotingHouseholderQR& compute(const MatrixType& matrix); const IntRowVectorType& colsPermutation() const { - ei_assert(m_isInitialized && "RRQR is not initialized."); + ei_assert(m_isInitialized && "FULLPIVOTINGHOUSEHOLDERQR is not initialized."); return m_cols_permutation; } const IntColVectorType& rowsTranspositions() const { - ei_assert(m_isInitialized && "RRQR is not initialized."); + ei_assert(m_isInitialized && "FULLPIVOTINGHOUSEHOLDERQR is not initialized."); return m_rows_transpositions; } inline int rank() const { - ei_assert(m_isInitialized && "RRQR is not initialized."); + ei_assert(m_isInitialized && "FULLPIVOTINGHOUSEHOLDERQR is not initialized."); return m_rank; } @@ -136,7 +136,7 @@ template class HouseholderRRQR #ifndef EIGEN_HIDE_HEAVY_CODE template -HouseholderRRQR& HouseholderRRQR::compute(const MatrixType& matrix) +FullPivotingHouseholderQR& FullPivotingHouseholderQR::compute(const MatrixType& matrix) { int rows = matrix.rows(); int cols = matrix.cols(); @@ -148,7 +148,6 @@ HouseholderRRQR& HouseholderRRQR::compute(const MatrixTy RowVectorType temp(cols); - // TODO: experiment to see the best formula m_precision = epsilon() * size; m_rows_transpositions.resize(matrix.rows()); @@ -198,7 +197,6 @@ HouseholderRRQR& HouseholderRRQR::compute(const MatrixTy m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); m_qr.coeffRef(k,k) = beta; - // apply H to remaining part of m_qr from the left m_qr.corner(BottomRight, rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); } @@ -215,12 +213,12 @@ HouseholderRRQR& HouseholderRRQR::compute(const MatrixTy template template -void HouseholderRRQR::solve( +void FullPivotingHouseholderQR::solve( const MatrixBase& b, ResultType *result ) const { - ei_assert(m_isInitialized && "HouseholderRRQR is not initialized."); + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); @@ -247,9 +245,9 @@ void HouseholderRRQR::solve( /** \returns the matrix Q */ template -MatrixType HouseholderRRQR::matrixQ() const +MatrixType FullPivotingHouseholderQR::matrixQ() const { - ei_assert(m_isInitialized && "HouseholderRRQR is not initialized."); + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); // compute the product H'_0 H'_1 ... H'_n-1, // where H_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), ...] @@ -269,18 +267,15 @@ MatrixType HouseholderRRQR::matrixQ() const #endif // EIGEN_HIDE_HEAVY_CODE -#if 0 -/** \return the Householder QR decomposition of \c *this. +/** \return the full-pivoting Householder QR decomposition of \c *this. * - * \sa class HouseholderRRQR + * \sa class FullPivotingHouseholderQR */ template -const HouseholderRRQR::PlainMatrixType> -MatrixBase::householderQr() const +const FullPivotingHouseholderQR::PlainMatrixType> +MatrixBase::fullPivotingHouseholderQr() const { - return HouseholderRRQR(eval()); + return FullPivotingHouseholderQR(eval()); } -#endif - -#endif // EIGEN_QR_H +#endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b79e069b7..fc43bbb1d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -121,7 +121,8 @@ ei_add_test(lu ${EI_OFLAG}) ei_add_test(determinant) ei_add_test(inverse ${EI_OFLAG}) ei_add_test(qr) -ei_add_test(rrqr) +ei_add_test(qr_colpivoting) +ei_add_test(qr_fullpivoting) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") ei_add_test(svd) diff --git a/test/rrqr.cpp b/test/qr_colpivoting.cpp similarity index 95% rename from test/rrqr.cpp rename to test/qr_colpivoting.cpp index b6cc75d17..3a07c7131 100644 --- a/test/rrqr.cpp +++ b/test/qr_colpivoting.cpp @@ -37,7 +37,7 @@ template void qr() typedef Matrix VectorType; MatrixType m1; createRandomMatrixOfRank(rank,rows,cols,m1); - HouseholderRRQR qr(m1); + ColPivotingHouseholderQR qr(m1); VERIFY_IS_APPROX(rank, qr.rank()); MatrixType r = qr.matrixQR(); @@ -74,7 +74,7 @@ template void qr_invertible() m1 += a * a.adjoint(); } - HouseholderRRQR qr(m1); + ColPivotingHouseholderQR qr(m1); m3 = MatrixType::Random(size,size); qr.solve(m3, &m2); VERIFY_IS_APPROX(m3, m1*m2); @@ -84,13 +84,13 @@ template void qr_verify_assert() { MatrixType tmp; - HouseholderRRQR qr; + ColPivotingHouseholderQR qr; VERIFY_RAISES_ASSERT(qr.matrixR()) VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) VERIFY_RAISES_ASSERT(qr.matrixQ()) } -void test_rrqr() +void test_qr_colpivoting() { for(int i = 0; i < 1; i++) { // FIXME : very weird bug here diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp new file mode 100644 index 000000000..fcbf08e4c --- /dev/null +++ b/test/qr_fullpivoting.cpp @@ -0,0 +1,115 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2009 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 . + +#include "main.h" +#include + +template void qr() +{ + /* this test covers the following files: QR.h */ + int rows = ei_random(20,200), cols = ei_random(20,200), cols2 = ei_random(20,200); + int rank = ei_random(1, std::min(rows, cols)-1); + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix SquareMatrixType; + typedef Matrix VectorType; + MatrixType m1; + createRandomMatrixOfRank(rank,rows,cols,m1); + FullPivotingHouseholderQR qr(m1); + VERIFY_IS_APPROX(rank, qr.rank()); + + MatrixType r = qr.matrixQR(); + // FIXME need better way to construct trapezoid + for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); + + MatrixType b = qr.matrixQ() * r; + + MatrixType c = MatrixType::Zero(rows,cols); + + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); + VERIFY_IS_APPROX(m1, c); + + MatrixType m2 = MatrixType::Random(cols,cols2); + MatrixType m3 = m1*m2; + m2 = MatrixType::Random(cols,cols2); + qr.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); +} + +template void qr_invertible() +{ + typedef typename NumTraits::Real RealScalar; + int size = ei_random(10,50); + + MatrixType m1(size, size), m2(size, size), m3(size, size); + m1 = MatrixType::Random(size,size); + + if (ei_is_same_type::ret) + { + // let's build a matrix more stable to inverse + MatrixType a = MatrixType::Random(size,size*2); + m1 += a * a.adjoint(); + } + + FullPivotingHouseholderQR qr(m1); + m3 = MatrixType::Random(size,size); + qr.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); +} + +template void qr_verify_assert() +{ + MatrixType tmp; + + FullPivotingHouseholderQR qr; + VERIFY_RAISES_ASSERT(qr.matrixR()) + VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(qr.matrixQ()) +} + +void test_qr_fullpivoting() +{ + for(int i = 0; i < 1; i++) { + // FIXME : very weird bug here +// CALL_SUBTEST( qr(Matrix2f()) ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + CALL_SUBTEST( qr() ); + } + + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + CALL_SUBTEST( qr_invertible() ); + } + + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); +} From 09265496595c0b4abf2260332b26f086caf8efa1 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 24 Aug 2009 00:02:22 -0400 Subject: [PATCH 12/16] fix bug: with complex matrices, the condition (ei_imag(c0)==RealScalar(0)) being wrong could bypass the other condition in the &&. at least that's my explanation why the test_lu was often failing on complex matrices (it uses that via createRandomMatrixOfRank) and why that's fixed by this diff. also gcc 4.4 gave a warning about tailSqNorm potentially uninitialized --- Eigen/src/Householder/Householder.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 8a274d240..ea39e4c30 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -74,10 +74,10 @@ void MatrixBase::makeHouseholder( EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) VectorBlock tail(derived(), 1, size()-1); - RealScalar tailSqNorm; + RealScalar tailSqNorm = size()==1 ? 0 : tail.squaredNorm(); Scalar c0 = coeff(0); - if( (size()==1 || (tailSqNorm=tail.squaredNorm()) == RealScalar(0)) && ei_imag(c0)==RealScalar(0)) + if( tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0)) { *tau = 0; *beta = ei_real(c0); From b37ab9b324ca55e365146266fc93372910fb0a98 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 24 Aug 2009 00:02:49 -0400 Subject: [PATCH 13/16] clarifications in LU::solve() and in LU documentation --- Eigen/src/LU/LU.h | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index d32fe08a1..d4f5fce1a 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -42,11 +42,10 @@ * This decomposition provides the generic approach to solving systems of linear equations, computing * the rank, invertibility, inverse, kernel, and determinant. * - * This LU decomposition is very stable and well tested with large matrices. Even exact rank computation - * works at sizes larger than 1000x1000. However there are use cases where the SVD decomposition is inherently - * more stable when dealing with numerically damaged input. For example, computing the kernel is more stable with - * SVD because the SVD can determine which singular values are negligible while LU has to work at the level of matrix - * coefficients that are less meaningful in this respect. + * This LU decomposition is very stable and well tested with large matrices. However there are use cases where the SVD + * decomposition is inherently more stable and/or flexible. For example, when computing the kernel of a matrix, + * working with the SVD allows to select the smallest singular values of the matrix, something that + * the LU decomposition doesn't see. * * The data of the LU decomposition can be directly accessed through the methods matrixLU(), * permutationP(), permutationQ(). @@ -326,6 +325,7 @@ template class LU inline void computeInverse(MatrixType *result) const { ei_assert(m_originalMatrix != 0 && "LU is not initialized."); + ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result); } @@ -456,6 +456,7 @@ template typename ei_traits::Scalar LU::determinant() const { ei_assert(m_originalMatrix != 0 && "LU is not initialized."); + ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!"); return Scalar(m_det_pq) * m_lu.diagonal().prod(); } @@ -533,7 +534,8 @@ bool LU::solve( ) const { ei_assert(m_originalMatrix != 0 && "LU is not initialized."); - + if(m_rank==0) return false; + /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. * So we proceed as follows: * Step 1: compute c = Pb. @@ -552,7 +554,8 @@ bool LU::solve( for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i); // Step 2 - m_lu.corner(Eigen::TopLeft,smalldim,smalldim).template triangularView() + m_lu.corner(Eigen::TopLeft,smalldim,smalldim) + .template triangularView() .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols())); if(rows>cols) { @@ -564,11 +567,10 @@ bool LU::solve( if(!isSurjective()) { // is c is in the image of U ? - RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : 0; - for(int col = 0; col < c.cols(); ++col) - for(int row = m_rank; row < c.rows(); ++row) - if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision)) - return false; + RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff(); + if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision)) + return false; } m_lu.corner(TopLeft, m_rank, m_rank) .template triangularView() From 154bdac9f4fe8b3b2b9a3089b87392c16ccd89e1 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 24 Aug 2009 00:09:01 -0400 Subject: [PATCH 14/16] small improvements --- test/lu.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/lu.cpp b/test/lu.cpp index 4ad92bb11..75680b96b 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -44,7 +44,7 @@ template void lu_non_invertible() VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); VERIFY(!lu.isInjective()); VERIFY(!lu.isInvertible()); - VERIFY(lu.isSurjective() == (lu.rank() == rows)); + VERIFY(!lu.isSurjective()); VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); VERIFY(m1image.lu().rank() == rank); MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); @@ -53,7 +53,7 @@ template void lu_non_invertible() m2 = MatrixType::Random(cols,cols2); m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); - lu.solve(m3, &m2); + VERIFY(lu.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); m3 = MatrixType::Random(rows,cols2); VERIFY(!lu.solve(m3, &m2)); From c9a307f3303ba7b49bc119ee0043df88d53cd995 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 24 Aug 2009 00:23:35 -0400 Subject: [PATCH 15/16] give FullPivotingHouseholderQR all the modern comfort --- Eigen/src/Householder/Householder.h | 2 +- Eigen/src/QR/FullPivotingHouseholderQR.h | 128 +++++++++++++++++++++-- test/qr_fullpivoting.cpp | 15 ++- 3 files changed, 135 insertions(+), 10 deletions(-) diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index ea39e4c30..36f02d7ce 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -77,7 +77,7 @@ void MatrixBase::makeHouseholder( RealScalar tailSqNorm = size()==1 ? 0 : tail.squaredNorm(); Scalar c0 = coeff(0); - if( tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0)) + if(tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0)) { *tau = 0; *beta = ei_real(c0); diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index f7b0f1cc1..0ffcfe88c 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -94,7 +94,7 @@ template class FullPivotingHouseholderQR * Output: \verbinclude FullPivotingHouseholderQR_solve.out */ template - void solve(const MatrixBase& b, ResultType *result) const; + bool solve(const MatrixBase& b, ResultType *result) const; MatrixType matrixQ(void) const; @@ -106,22 +106,117 @@ template class FullPivotingHouseholderQR const IntRowVectorType& colsPermutation() const { - ei_assert(m_isInitialized && "FULLPIVOTINGHOUSEHOLDERQR is not initialized."); + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); return m_cols_permutation; } const IntColVectorType& rowsTranspositions() const { - ei_assert(m_isInitialized && "FULLPIVOTINGHOUSEHOLDERQR is not initialized."); + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); return m_rows_transpositions; } + /** \returns the absolute value of the determinant of the matrix of which + * *this is the QR decomposition. It has only linear complexity + * (that is, O(n) where n is the dimension of the square matrix) + * as the QR decomposition has already been computed. + * + * \note This is only for square matrices. + * + * \warning a determinant can be very big or small, so for matrices + * of large enough dimension, there is a risk of overflow/underflow. + * + * \sa MatrixBase::determinant() + */ + typename MatrixType::RealScalar absDeterminant() const; + + /** \returns the rank of the matrix of which *this is the QR decomposition. + * + * \note This is computed at the time of the construction of the QR decomposition. This + * method does not perform any further computation. + */ inline int rank() const { - ei_assert(m_isInitialized && "FULLPIVOTINGHOUSEHOLDERQR is not initialized."); + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); return m_rank; } + /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. + * + * \note Since the rank is computed at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline int dimensionOfKernel() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_qr.cols() - m_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 at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline bool isInjective() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_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 at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline bool isSurjective() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return m_rank == m_qr.rows(); + } + + /** \returns true if the matrix of which *this is the QR decomposition is invertible. + * + * \note Since the rank is computed at the time of the construction of the QR decomposition, this + * method almost does not perform any further computation. + */ + inline bool isInvertible() const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + return isInjective() && isSurjective(); + } + + /** Computes the inverse of the matrix of which *this is the QR decomposition. + * + * \param result a pointer to the matrix into which to store the inverse. Resized if needed. + * + * \note If this matrix is not invertible, *result is left with undefined coefficients. + * Use isInvertible() to first determine whether this matrix is invertible. + * + * \sa inverse() + */ + inline void computeInverse(MatrixType *result) const + { + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!"); + solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result); + } + + /** \returns the inverse of the matrix of which *this is the QR decomposition. + * + * \note If this matrix is not invertible, the returned matrix has undefined coefficients. + * Use isInvertible() to first determine whether this matrix is invertible. + * + * \sa computeInverse() + */ + inline MatrixType inverse() const + { + MatrixType result; + computeInverse(&result); + return result; + } + protected: MatrixType m_qr; HCoeffsType m_hCoeffs; @@ -135,6 +230,14 @@ template class FullPivotingHouseholderQR #ifndef EIGEN_HIDE_HEAVY_CODE +template +typename MatrixType::RealScalar FullPivotingHouseholderQR::absDeterminant() const +{ + ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); + return ei_abs(m_qr.diagonal().prod()); +} + template FullPivotingHouseholderQR& FullPivotingHouseholderQR::compute(const MatrixType& matrix) { @@ -163,8 +266,8 @@ FullPivotingHouseholderQR& FullPivotingHouseholderQR::co RealScalar biggest_in_corner; biggest_in_corner = m_qr.corner(Eigen::BottomRight, rows-k, cols-k) - .cwise().abs() - .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); + .cwise().abs() + .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); row_of_biggest_in_corner += k; col_of_biggest_in_corner += k; if(k==0) biggest = biggest_in_corner; @@ -213,12 +316,14 @@ FullPivotingHouseholderQR& FullPivotingHouseholderQR::co template template -void FullPivotingHouseholderQR::solve( +bool FullPivotingHouseholderQR::solve( const MatrixBase& b, ResultType *result ) const { ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); + if(m_rank==0) return false; + const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); @@ -234,6 +339,14 @@ void FullPivotingHouseholderQR::solve( .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); } + if(!isSurjective()) + { + // is c is in the image of R ? + RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff(); + if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision)) + return false; + } m_qr.corner(TopLeft, m_rank, m_rank) .template triangularView() .solveInPlace(c.corner(TopLeft, m_rank, c.cols())); @@ -241,6 +354,7 @@ void FullPivotingHouseholderQR::solve( result->resize(m_qr.cols(), b.cols()); for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i); for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero(); + return true; } /** \returns the matrix Q */ diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index fcbf08e4c..1e0601a83 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -39,6 +39,11 @@ template void qr() createRandomMatrixOfRank(rank,rows,cols,m1); FullPivotingHouseholderQR qr(m1); VERIFY_IS_APPROX(rank, qr.rank()); + VERIFY(cols - qr.rank() == qr.dimensionOfKernel()); + VERIFY(!qr.isInjective()); + VERIFY(!qr.isInvertible()); + VERIFY(!qr.isSurjective()); + MatrixType r = qr.matrixQR(); // FIXME need better way to construct trapezoid @@ -54,8 +59,10 @@ template void qr() MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); - qr.solve(m3, &m2); + VERIFY(qr.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); + m3 = MatrixType::Random(rows,cols2); + VERIFY(!qr.solve(m3, &m2)); } template void qr_invertible() @@ -74,8 +81,12 @@ template void qr_invertible() } FullPivotingHouseholderQR qr(m1); + VERIFY(qr.isInjective()); + VERIFY(qr.isInvertible()); + VERIFY(qr.isSurjective()); + m3 = MatrixType::Random(size,size); - qr.solve(m3, &m2); + VERIFY(qr.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); } From f31b5a71148387310a55a96158a494e83a19a0e2 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 24 Aug 2009 00:35:42 -0400 Subject: [PATCH 16/16] add test for absDeterminant() --- test/qr_fullpivoting.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index 1e0601a83..d784e0d43 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -68,6 +68,8 @@ template void qr() template void qr_invertible() { typedef typename NumTraits::Real RealScalar; + typedef typename MatrixType::Scalar Scalar; + int size = ei_random(10,50); MatrixType m1(size, size), m2(size, size), m3(size, size); @@ -88,6 +90,15 @@ template void qr_invertible() m3 = MatrixType::Random(size,size); VERIFY(qr.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); + + // now construct a matrix with prescribed determinant + m1.setZero(); + for(int i = 0; i < size; i++) m1(i,i) = ei_random(); + RealScalar absdet = ei_abs(m1.diagonal().prod()); + m3 = qr.matrixQ(); // get a unitary + m1 = m3 * m1 * m3; + qr.compute(m1); + VERIFY_IS_APPROX(absdet, qr.absDeterminant()); } template void qr_verify_assert()