diff options
| -rw-r--r-- | Eigen/CMakeLists.txt | 20 | ||||
| -rw-r--r-- | Eigen/LU | 3 | ||||
| -rw-r--r-- | Eigen/src/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/ReturnByValue.h | 6 | ||||
| -rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 7 | ||||
| -rw-r--r-- | Eigen/src/LU/FullPivLU.h | 257 | ||||
| -rw-r--r-- | Eigen/src/misc/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | Eigen/src/misc/Image.h | 72 | ||||
| -rw-r--r-- | Eigen/src/misc/Kernel.h | 71 | ||||
| -rw-r--r-- | Eigen/src/misc/Solve.h | 65 | ||||
| -rw-r--r-- | test/lu.cpp | 8 |
11 files changed, 329 insertions, 187 deletions
diff --git a/Eigen/CMakeLists.txt b/Eigen/CMakeLists.txt index 931cc6e20..e0eb837a5 100644 --- a/Eigen/CMakeLists.txt +++ b/Eigen/CMakeLists.txt @@ -1,25 +1,5 @@ set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi Eigenvalues) -if(EIGEN_BUILD_LIB) - set(Eigen_SRCS - src/Core/CoreInstantiations.cpp - src/Cholesky/CholeskyInstantiations.cpp - src/QR/QrInstantiations.cpp - ) - - add_library(Eigen2 SHARED ${Eigen_SRCS}) - - install(TARGETS Eigen2 - RUNTIME DESTINATION bin - LIBRARY DESTINATION lib - ARCHIVE DESTINATION lib) -endif(EIGEN_BUILD_LIB) - -if(CMAKE_COMPILER_IS_GNUCXX) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g1 -O2") - set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -g1 -O2") -endif(CMAKE_COMPILER_IS_GNUCXX) - install(FILES ${Eigen_HEADERS} DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen COMPONENT Devel @@ -18,6 +18,9 @@ namespace Eigen { * \endcode */ +#include "src/misc/Solve.h" +#include "src/misc/Kernel.h" +#include "src/misc/Image.h" #include "src/LU/FullPivLU.h" #include "src/LU/PartialPivLU.h" #include "src/LU/Determinant.h" diff --git a/Eigen/src/CMakeLists.txt b/Eigen/src/CMakeLists.txt index 0df8273d1..4c3ee6cb8 100644 --- a/Eigen/src/CMakeLists.txt +++ b/Eigen/src/CMakeLists.txt @@ -10,3 +10,4 @@ ADD_SUBDIRECTORY(Sparse) ADD_SUBDIRECTORY(Jacobi) ADD_SUBDIRECTORY(Householder) ADD_SUBDIRECTORY(Eigenvalues) +ADD_SUBDIRECTORY(misc)
\ No newline at end of file diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h index 87b057f86..1d977ace2 100644 --- a/Eigen/src/Core/ReturnByValue.h +++ b/Eigen/src/Core/ReturnByValue.h @@ -38,8 +38,10 @@ struct ei_traits<ReturnByValue<Derived> > // matrix.inverse().block(...) // because the Block ctor with direct access // wants to call coeffRef() to get an address, and that fails (infinite recursion) as ReturnByValue - // doesnt implement coeffRef(). The better fix is probably rather to make Block work directly - // on the nested type, right? + // doesnt implement coeffRef(). + // The fact that I had to do that shows that when doing xpr.block() with a non-direct-access xpr, + // even if xpr has the EvalBeforeNestingBit, the block() doesn't use direct access on the evaluated + // xpr. Flags = (ei_traits<typename ei_traits<Derived>::ReturnMatrixType>::Flags | EvalBeforeNestingBit) & ~DirectAccessBit }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index c8f2c4cd7..b54f71ed1 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -66,6 +66,13 @@ template<typename ExpressionType> class WithFormat; template<typename MatrixType> struct CommaInitializer; template<typename Derived> class ReturnByValue; +template<typename DecompositionType, typename Rhs> struct ei_solve_return_value; +template<typename DecompositionType, typename Rhs, typename Dest> struct ei_solve_impl; +template<typename DecompositionType> struct ei_kernel_return_value; +template<typename DecompositionType, typename Dest> struct ei_kernel_impl; +template<typename DecompositionType> struct ei_image_return_value; +template<typename DecompositionType, typename Dest> struct ei_image_impl; + template<typename _Scalar, int Rows=Dynamic, int Cols=Dynamic, int Supers=Dynamic, int Subs=Dynamic, int Options=0> class BandMatrix; diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 067b59549..c5f44dcea 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -25,10 +25,6 @@ #ifndef EIGEN_LU_H #define EIGEN_LU_H -template<typename MatrixType, typename Rhs> struct ei_fullpivlu_solve_impl; -template<typename MatrixType> struct ei_fullpivlu_kernel_impl; -template<typename MatrixType> struct ei_fullpivlu_image_impl; - /** \ingroup LU_Module * * \class FullPivLU @@ -59,10 +55,10 @@ template<typename MatrixType> struct ei_fullpivlu_image_impl; * * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse() */ -template<typename MatrixType> class FullPivLU +template<typename _MatrixType> class FullPivLU { public: - + typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; @@ -70,17 +66,12 @@ template<typename MatrixType> class FullPivLU typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType; - enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( - MatrixType::MaxColsAtCompileTime, - MatrixType::MaxRowsAtCompileTime) - }; - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via LU::compute(const MatrixType&). - */ + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via LU::compute(const MatrixType&). + */ FullPivLU(); /** Constructor. @@ -167,10 +158,10 @@ template<typename MatrixType> class FullPivLU * * \sa image() */ - inline const ei_fullpivlu_kernel_impl<MatrixType> kernel() const + inline const ei_kernel_return_value<FullPivLU> kernel() const { ei_assert(m_isInitialized && "LU is not initialized."); - return ei_fullpivlu_kernel_impl<MatrixType>(*this); + return ei_kernel_return_value<FullPivLU>(*this); } /** \returns the image of the matrix, also called its column-space. The columns of the returned matrix @@ -192,12 +183,11 @@ template<typename MatrixType> class FullPivLU * * \sa kernel() */ - template<typename OriginalMatrixType> - inline const ei_fullpivlu_image_impl<MatrixType> - image(const MatrixBase<OriginalMatrixType>& originalMatrix) const + inline const ei_image_return_value<FullPivLU> + image(const MatrixType& originalMatrix) const { ei_assert(m_isInitialized && "LU is not initialized."); - return ei_fullpivlu_image_impl<MatrixType>(*this, originalMatrix.derived()); + return ei_image_return_value<FullPivLU>(*this, originalMatrix); } /** \return a solution x to the equation Ax=b, where A is the matrix of which @@ -220,11 +210,11 @@ template<typename MatrixType> class FullPivLU * \sa TriangularView::solve(), kernel(), inverse() */ template<typename Rhs> - inline const ei_fullpivlu_solve_impl<MatrixType, Rhs> + inline const ei_solve_return_value<FullPivLU, Rhs> solve(const MatrixBase<Rhs>& b) const { ei_assert(m_isInitialized && "LU is not initialized."); - return ei_fullpivlu_solve_impl<MatrixType, Rhs>(*this, b.derived()); + return ei_solve_return_value<FullPivLU, Rhs>(*this, b.derived()); } /** \returns the determinant of the matrix of which @@ -365,14 +355,17 @@ template<typename MatrixType> class FullPivLU * * \sa MatrixBase::inverse() */ - inline const ei_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const + inline const ei_solve_return_value<FullPivLU,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const { ei_assert(m_isInitialized && "LU is not initialized."); ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); - return ei_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > + return ei_solve_return_value<FullPivLU,NestByValue<typename MatrixType::IdentityReturnType> > (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); } + inline int rows() const { return m_lu.rows(); } + inline int cols() const { return m_lu.cols(); } + protected: MatrixType m_lu; IntColVectorType m_p; @@ -492,42 +485,21 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons /********* Implementation of kernel() **************************************************/ -template<typename MatrixType> -struct ei_traits<ei_fullpivlu_kernel_impl<MatrixType> > -{ - typedef Matrix< - typename MatrixType::Scalar, - MatrixType::ColsAtCompileTime, // the number of rows in the "kernel matrix" - // is the number of cols of the original matrix - // so that the product "matrix * kernel = zero" makes sense - Dynamic, // we don't know at compile-time the dimension of the kernel - MatrixType::Options, - MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter - MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space, - // whose dimension is the number of columns of the original matrix - > ReturnMatrixType; -}; - -template<typename MatrixType> -struct ei_fullpivlu_kernel_impl : public ReturnByValue<ei_fullpivlu_kernel_impl<MatrixType> > +template<typename MatrixType, typename Dest> +struct ei_kernel_impl<FullPivLU<MatrixType>, Dest> + : ei_kernel_return_value<FullPivLU<MatrixType> > { - typedef FullPivLU<MatrixType> LUType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - const LUType& m_lu; - int m_rank, m_cols; - - ei_fullpivlu_kernel_impl(const LUType& lu) - : m_lu(lu), - m_rank(lu.rank()), - m_cols(m_rank==lu.matrixLU().cols() ? 1 : lu.matrixLU().cols() - m_rank){} - - inline int rows() const { return m_lu.matrixLU().cols(); } - inline int cols() const { return m_cols; } - - template<typename Dest> void evalTo(Dest& dst) const + enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( + MatrixType::MaxColsAtCompileTime, + MatrixType::MaxRowsAtCompileTime) + }; + + void evalTo(Dest& dst) const { - const int cols = m_lu.matrixLU().cols(), dimker = cols - m_rank; + const FullPivLU<MatrixType>& dec = this->m_dec; + const int cols = dec.matrixLU().cols(), rank = this->m_rank, dimker = cols - rank; if(dimker == 0) { // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's @@ -549,89 +521,73 @@ struct ei_fullpivlu_kernel_impl : public ReturnByValue<ei_fullpivlu_kernel_impl< * * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end. * Thus, the diagonal of U ends with exactly - * m_dimKer zero's. Let us use that to construct dimKer linearly + * dimKer zero's. Let us use that to construct dimKer linearly * independent vectors in Ker U. */ - Matrix<int, Dynamic, 1, 0, LUType::MaxSmallDimAtCompileTime, 1> pivots(m_rank); - RealScalar premultiplied_threshold = m_lu.maxPivot() * m_lu.threshold(); + Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank); + RealScalar premultiplied_threshold = dec.maxPivot() * dec.threshold(); int p = 0; - for(int i = 0; i < m_lu.nonzeroPivots(); ++i) - if(ei_abs(m_lu.matrixLU().coeff(i,i)) > premultiplied_threshold) + for(int i = 0; i < dec.nonzeroPivots(); ++i) + if(ei_abs(dec.matrixLU().coeff(i,i)) > premultiplied_threshold) pivots.coeffRef(p++) = i; - ei_assert(p == m_rank && "You hit a bug in Eigen! Please report (backtrace and matrix)!"); + ei_internal_assert(p == rank); // we construct a temporaty trapezoid matrix m, by taking the U matrix and // permuting the rows and cols to bring the nonnegligible pivots to the top of // the main diagonal. We need that to be able to apply our triangular solvers. // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options, - LUType::MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime> - m(m_lu.matrixLU().block(0, 0, m_rank, cols)); - for(int i = 0; i < m_rank; ++i) + MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime> + m(dec.matrixLU().block(0, 0, rank, cols)); + for(int i = 0; i < rank; ++i) { if(i) m.row(i).start(i).setZero(); - m.row(i).end(cols-i) = m_lu.matrixLU().row(pivots.coeff(i)).end(cols-i); + m.row(i).end(cols-i) = dec.matrixLU().row(pivots.coeff(i)).end(cols-i); } - m.block(0, 0, m_rank, m_rank).template triangularView<StrictlyLowerTriangular>().setZero(); - for(int i = 0; i < m_rank; ++i) + m.block(0, 0, rank, rank); + m.block(0, 0, rank, rank).template triangularView<StrictlyLowerTriangular>().setZero(); + for(int i = 0; i < rank; ++i) m.col(i).swap(m.col(pivots.coeff(i))); // ok, we have our trapezoid matrix, we can apply the triangular solver. // notice that the math behind this suggests that we should apply this to the // negative of the RHS, but for performance we just put the negative sign elsewhere, see below. - m.corner(TopLeft, m_rank, m_rank) + m.corner(TopLeft, rank, rank) .template triangularView<UpperTriangular>().solveInPlace( - m.corner(TopRight, m_rank, dimker) + m.corner(TopRight, rank, dimker) ); // now we must undo the column permutation that we had applied! - for(int i = m_rank-1; i >= 0; --i) + for(int i = rank-1; i >= 0; --i) m.col(i).swap(m.col(pivots.coeff(i))); // see the negative sign in the next line, that's what we were talking about above. - for(int i = 0; i < m_rank; ++i) dst.row(m_lu.permutationQ().coeff(i)) = -m.row(i).end(dimker); - for(int i = m_rank; i < cols; ++i) dst.row(m_lu.permutationQ().coeff(i)).setZero(); - for(int k = 0; k < dimker; ++k) dst.coeffRef(m_lu.permutationQ().coeff(m_rank+k), k) = Scalar(1); + for(int i = 0; i < rank; ++i) dst.row(dec.permutationQ().coeff(i)) = -m.row(i).end(dimker); + for(int i = rank; i < cols; ++i) dst.row(dec.permutationQ().coeff(i)).setZero(); + for(int k = 0; k < dimker; ++k) dst.coeffRef(dec.permutationQ().coeff(rank+k), k) = Scalar(1); } }; /***** Implementation of image() *****************************************************/ -template<typename MatrixType> -struct ei_traits<ei_fullpivlu_image_impl<MatrixType> > -{ - typedef Matrix< - typename MatrixType::Scalar, - MatrixType::RowsAtCompileTime, // the image is a subspace of the destination space, whose - // dimension is the number of rows of the original matrix - Dynamic, // we don't know at compile time the dimension of the image (the rank) - MatrixType::Options, - MatrixType::MaxRowsAtCompileTime, // the image matrix will consist of columns from the original matrix, - MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns. - > ReturnMatrixType; -}; - -template<typename MatrixType> -struct ei_fullpivlu_image_impl : public ReturnByValue<ei_fullpivlu_image_impl<MatrixType> > +template<typename MatrixType, typename Dest> +struct ei_image_impl<FullPivLU<MatrixType>, Dest> + : ei_image_return_value<FullPivLU<MatrixType> > { - typedef FullPivLU<MatrixType> LUType; + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - const LUType& m_lu; - int m_rank, m_cols; - const MatrixType& m_originalMatrix; - - ei_fullpivlu_image_impl(const LUType& lu, const MatrixType& originalMatrix) - : m_lu(lu), m_rank(lu.rank()), - m_cols(m_rank == 0 ? 1 : m_rank), - m_originalMatrix(originalMatrix) {} - - inline int rows() const { return m_lu.matrixLU().rows(); } - inline int cols() const { return m_cols; } - - template<typename Dest> void evalTo(Dest& dst) const + enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( + MatrixType::MaxColsAtCompileTime, + MatrixType::MaxRowsAtCompileTime) + }; + + void evalTo(Dest& dst) const { - if(m_rank == 0) + const int rank = this->m_rank; + const FullPivLU<MatrixType>& dec = this->m_dec; + const MatrixType& originalMatrix = this->m_originalMatrix; + if(rank == 0) { // The Image is just {0}, so it doesn't have a basis properly speaking, but let's // avoid crashing/asserting as that depends on floating point calculations. Let's @@ -640,61 +596,40 @@ struct ei_fullpivlu_image_impl : public ReturnByValue<ei_fullpivlu_image_impl<Ma return; } - Matrix<int, Dynamic, 1, 0, LUType::MaxSmallDimAtCompileTime, 1> pivots(m_rank); - RealScalar premultiplied_threshold = m_lu.maxPivot() * m_lu.threshold(); + Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank); + RealScalar premultiplied_threshold = dec.maxPivot() * dec.threshold(); int p = 0; - for(int i = 0; i < m_lu.nonzeroPivots(); ++i) - if(ei_abs(m_lu.matrixLU().coeff(i,i)) > premultiplied_threshold) + for(int i = 0; i < dec.nonzeroPivots(); ++i) + if(ei_abs(dec.matrixLU().coeff(i,i)) > premultiplied_threshold) pivots.coeffRef(p++) = i; - ei_assert(p == m_rank && "You hit a bug in Eigen! Please report (backtrace and matrix)!"); + ei_internal_assert(p == rank); - for(int i = 0; i < m_rank; ++i) - dst.col(i) = m_originalMatrix.col(m_lu.permutationQ().coeff(pivots.coeff(i))); + for(int i = 0; i < rank; ++i) + dst.col(i) = originalMatrix.col(dec.permutationQ().coeff(pivots.coeff(i))); } }; /***** Implementation of solve() *****************************************************/ -template<typename MatrixType,typename Rhs> -struct ei_traits<ei_fullpivlu_solve_impl<MatrixType,Rhs> > -{ - typedef Matrix<typename Rhs::Scalar, - MatrixType::ColsAtCompileTime, - Rhs::ColsAtCompileTime, - Rhs::PlainMatrixType::Options, - MatrixType::MaxColsAtCompileTime, - Rhs::MaxColsAtCompileTime> ReturnMatrixType; -}; - -template<typename MatrixType, typename Rhs> -struct ei_fullpivlu_solve_impl : public ReturnByValue<ei_fullpivlu_solve_impl<MatrixType, Rhs> > +template<typename MatrixType, typename Rhs, typename Dest> +struct ei_solve_impl<FullPivLU<MatrixType>, Rhs, Dest> + : ei_solve_return_value<FullPivLU<MatrixType>, Rhs> { - typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested; - typedef FullPivLU<MatrixType> LUType; - const LUType& m_lu; - const typename Rhs::Nested m_rhs; - - ei_fullpivlu_solve_impl(const LUType& lu, const Rhs& rhs) - : m_lu(lu), m_rhs(rhs) - {} - - inline int rows() const { return m_lu.matrixLU().cols(); } - inline int cols() const { return m_rhs.cols(); } - - template<typename Dest> void evalTo(Dest& dst) const + void evalTo(Dest& dst) const { /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. - * So we proceed as follows: - * Step 1: compute c = P * rhs. - * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible. - * Step 3: replace c by the solution x to Ux = c. May or may not exist. - * Step 4: result = Q * c; - */ - - const int rows = m_lu.matrixLU().rows(), - cols = m_lu.matrixLU().cols(), - nonzero_pivots = m_lu.nonzeroPivots(); - ei_assert(m_rhs.rows() == rows); + * So we proceed as follows: + * Step 1: compute c = P * rhs. + * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible. + * Step 3: replace c by the solution x to Ux = c. May or may not exist. + * Step 4: result = Q * c; + */ + + const FullPivLU<MatrixType>& dec = this->m_dec; + const Rhs& rhs = this->m_rhs; + const int rows = dec.matrixLU().rows(), cols = dec.matrixLU().cols(), + nonzero_pivots = dec.nonzeroPivots(); + ei_assert(rhs.rows() == rows); const int smalldim = std::min(rows, cols); if(nonzero_pivots == 0) @@ -702,36 +637,36 @@ struct ei_fullpivlu_solve_impl : public ReturnByValue<ei_fullpivlu_solve_impl<Ma dst.setZero(); return; } - - typename Rhs::PlainMatrixType c(m_rhs.rows(), m_rhs.cols()); + + typename Rhs::PlainMatrixType c(rhs.rows(), rhs.cols()); // Step 1 for(int i = 0; i < rows; ++i) - c.row(m_lu.permutationP().coeff(i)) = m_rhs.row(i); + c.row(dec.permutationP().coeff(i)) = rhs.row(i); // Step 2 - m_lu.matrixLU() + dec.matrixLU() .corner(Eigen::TopLeft,smalldim,smalldim) .template triangularView<UnitLowerTriangular>() .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols())); if(rows>cols) { c.corner(Eigen::BottomLeft, rows-cols, c.cols()) - -= m_lu.matrixLU().corner(Eigen::BottomLeft, rows-cols, cols) - * c.corner(Eigen::TopLeft, cols, c.cols()); + -= dec.matrixLU().corner(Eigen::BottomLeft, rows-cols, cols) + * c.corner(Eigen::TopLeft, cols, c.cols()); } // Step 3 - m_lu.matrixLU() + dec.matrixLU() .corner(TopLeft, nonzero_pivots, nonzero_pivots) .template triangularView<UpperTriangular>() .solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); // Step 4 for(int i = 0; i < nonzero_pivots; ++i) - dst.row(m_lu.permutationQ().coeff(i)) = c.row(i); - for(int i = nonzero_pivots; i < m_lu.matrixLU().cols(); ++i) - dst.row(m_lu.permutationQ().coeff(i)).setZero(); + dst.row(dec.permutationQ().coeff(i)) = c.row(i); + for(int i = nonzero_pivots; i < dec.matrixLU().cols(); ++i) + dst.row(dec.permutationQ().coeff(i)).setZero(); } }; diff --git a/Eigen/src/misc/CMakeLists.txt b/Eigen/src/misc/CMakeLists.txt new file mode 100644 index 000000000..a58ffb745 --- /dev/null +++ b/Eigen/src/misc/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_misc_SRCS "*.h") + +INSTALL(FILES + ${Eigen_misc_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/misc COMPONENT Devel + ) diff --git a/Eigen/src/misc/Image.h b/Eigen/src/misc/Image.h new file mode 100644 index 000000000..a7e2bceec --- /dev/null +++ b/Eigen/src/misc/Image.h @@ -0,0 +1,72 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_MISC_IMAGE_H +#define EIGEN_MISC_IMAGE_H + +/** \class ei_image_return_value + * + */ +template<typename DecompositionType> +struct ei_traits<ei_image_return_value<DecompositionType> > +{ + typedef typename DecompositionType::MatrixType MatrixType; + typedef Matrix< + typename MatrixType::Scalar, + MatrixType::RowsAtCompileTime, // the image is a subspace of the destination space, whose + // dimension is the number of rows of the original matrix + Dynamic, // we don't know at compile time the dimension of the image (the rank) + MatrixType::Options, + MatrixType::MaxRowsAtCompileTime, // the image matrix will consist of columns from the original matrix, + MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns. + > ReturnMatrixType; +}; + +template<typename _DecompositionType> struct ei_image_return_value + : public ReturnByValue<ei_image_return_value<_DecompositionType> > +{ + typedef _DecompositionType DecompositionType; + typedef typename DecompositionType::MatrixType MatrixType; + + const DecompositionType& m_dec; + int m_rank, m_cols; + const MatrixType& m_originalMatrix; + + ei_image_return_value(const DecompositionType& dec, const MatrixType& originalMatrix) + : m_dec(dec), m_rank(dec.rank()), + m_cols(m_rank == 0 ? 1 : m_rank), + m_originalMatrix(originalMatrix) + {} + + inline int rows() const { return m_dec.rows(); } + inline int cols() const { return m_cols; } + + template<typename Dest> inline void evalTo(Dest& dst) const + { + static_cast<const ei_image_impl<DecompositionType, Dest> *> + (this)->evalTo(dst); + } +}; + +#endif // EIGEN_MISC_IMAGE_H diff --git a/Eigen/src/misc/Kernel.h b/Eigen/src/misc/Kernel.h new file mode 100644 index 000000000..bfd75f54b --- /dev/null +++ b/Eigen/src/misc/Kernel.h @@ -0,0 +1,71 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_MISC_KERNEL_H +#define EIGEN_MISC_KERNEL_H + +/** \class ei_kernel_return_value + * + */ +template<typename DecompositionType> +struct ei_traits<ei_kernel_return_value<DecompositionType> > +{ + typedef typename DecompositionType::MatrixType MatrixType; + typedef Matrix< + typename MatrixType::Scalar, + MatrixType::ColsAtCompileTime, // the number of rows in the "kernel matrix" + // is the number of cols of the original matrix + // so that the product "matrix * kernel = zero" makes sense + Dynamic, // we don't know at compile-time the dimension of the kernel + MatrixType::Options, + MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter + MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space, + // whose dimension is the number of columns of the original matrix + > ReturnMatrixType; +}; + +template<typename _DecompositionType> struct ei_kernel_return_value + : public ReturnByValue<ei_kernel_return_value<_DecompositionType> > +{ + typedef _DecompositionType DecompositionType; + const DecompositionType& m_dec; + int m_rank, m_cols; + + ei_kernel_return_value(const DecompositionType& dec) + : m_dec(dec), + m_rank(dec.rank()), + m_cols(m_rank==dec.cols() ? 1 : dec.cols() - m_rank) + {} + + inline int rows() const { return m_dec.cols(); } + inline int cols() const { return m_cols; } + + template<typename Dest> inline void evalTo(Dest& dst) const + { + static_cast<const ei_kernel_impl<DecompositionType, Dest> *> + (this)->evalTo(dst); + } +}; + +#endif // EIGEN_MISC_KERNEL_H diff --git a/Eigen/src/misc/Solve.h b/Eigen/src/misc/Solve.h new file mode 100644 index 000000000..c62e34b13 --- /dev/null +++ b/Eigen/src/misc/Solve.h @@ -0,0 +1,65 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_MISC_SOLVE_H +#define EIGEN_MISC_SOLVE_H + +/** \class ei_solve_return_value + * + */ +template<typename DecompositionType, typename Rhs> +struct ei_traits<ei_solve_return_value<DecompositionType, Rhs> > +{ + typedef typename DecompositionType::MatrixType MatrixType; + typedef Matrix<typename Rhs::Scalar, + MatrixType::ColsAtCompileTime, + Rhs::ColsAtCompileTime, + Rhs::PlainMatrixType::Options, + MatrixType::MaxColsAtCompileTime, + Rhs::MaxColsAtCompileTime> ReturnMatrixType; +}; + +template<typename _DecompositionType, typename Rhs> struct ei_solve_return_value + : public ReturnByValue<ei_solve_return_value<_DecompositionType, Rhs> > +{ + typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNestedCleaned; + typedef _DecompositionType DecompositionType; + const DecompositionType& m_dec; + const typename Rhs::Nested m_rhs; + + ei_solve_return_value(const DecompositionType& dec, const Rhs& rhs) + : m_dec(dec), m_rhs(rhs) + {} + + inline int rows() const { return m_dec.cols(); } + inline int cols() const { return m_rhs.cols(); } + + template<typename Dest> inline void evalTo(Dest& dst) const + { + static_cast<const ei_solve_impl<DecompositionType, RhsNestedCleaned, Dest> *> + (this)->evalTo(dst); + } +}; + +#endif // EIGEN_MISC_SOLVE_H diff --git a/test/lu.cpp b/test/lu.cpp index c46ca9130..954893651 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -49,8 +49,8 @@ template<typename MatrixType> void lu_non_invertible() cols2 = cols = MatrixType::ColsAtCompileTime; } - typedef typename ei_fullpivlu_kernel_impl<MatrixType>::ReturnMatrixType KernelMatrixType; - typedef typename ei_fullpivlu_image_impl <MatrixType>::ReturnMatrixType ImageMatrixType; + typedef typename ei_kernel_return_value<FullPivLU<MatrixType> >::ReturnMatrixType KernelMatrixType; + typedef typename ei_image_return_value<FullPivLU<MatrixType> >::ReturnMatrixType ImageMatrixType; typedef Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> DynamicMatrixType; typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> CMatrixType; @@ -65,10 +65,10 @@ template<typename MatrixType> void lu_non_invertible() createRandomMatrixOfRank(rank, rows, cols, m1); FullPivLU<MatrixType> lu(m1); + std::cout << lu.kernel().rows() << " " << lu.kernel().cols() << std::endl; KernelMatrixType m1kernel = lu.kernel(); - ImageMatrixType m1image = lu.image(m1); + ImageMatrixType m1image = lu.image(m1); - // std::cerr << rank << " " << lu.rank() << std::endl; VERIFY(rank == lu.rank()); VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); VERIFY(!lu.isInjective()); |