From e4e58e8337e82ba76f6bf4fe7000acac9337056c Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 8 Nov 2009 16:51:41 -0500 Subject: [PATCH] simplifications in the ei_solve_impl system, factor out some boilerplate code --- Eigen/src/Cholesky/LDLT.h | 14 +-- Eigen/src/Cholesky/LLT.h | 15 +-- Eigen/src/Core/util/ForwardDeclarations.h | 6 +- Eigen/src/LU/FullPivLU.h | 114 +++++++++++----------- Eigen/src/LU/PartialPivLU.h | 25 +++-- Eigen/src/QR/ColPivHouseholderQR.h | 44 ++++----- Eigen/src/QR/FullPivHouseholderQR.h | 52 +++++----- Eigen/src/QR/HouseholderQR.h | 30 +++--- Eigen/src/SVD/SVD.h | 24 +++-- Eigen/src/misc/Image.h | 11 ++- Eigen/src/misc/Kernel.h | 10 +- Eigen/src/misc/Solve.h | 11 ++- doc/snippets/FullPivLU_image.cpp | 2 +- doc/snippets/FullPivLU_kernel.cpp | 2 +- doc/snippets/FullPivLU_solve.cpp | 2 +- doc/snippets/HouseholderQR_solve.cpp | 2 +- 16 files changed, 189 insertions(+), 175 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 761a4a8e8..a1faae49c 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -264,14 +264,16 @@ LDLT& LDLT::compute(const MatrixType& a) return *this; } -template -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const + EIGEN_MAKE_SOLVE_HELPERS(LDLT<_MatrixType>,Rhs) + + template void evalTo(Dest& dst) const { - dst = this->m_rhs; - this->m_dec.solveInPlace(dst); + dst = rhs(); + dec().solveInPlace(dst); } }; diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 0ad67aa5f..30c48578a 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -258,14 +258,17 @@ LLT& LLT::compute(const MatrixType& a) return *this; } -template -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const + typedef LLT<_MatrixType,UpLo> LLTType; + EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs) + + template void evalTo(Dest& dst) const { - dst = this->m_rhs; - this->m_dec.solveInPlace(dst); + dst = rhs(); + dec().solveInPlace(dst); } }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index b54f71ed1..323848919 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -67,11 +67,11 @@ template struct CommaInitializer; template class ReturnByValue; template struct ei_solve_return_value; -template struct ei_solve_impl; +template struct ei_solve_impl; template struct ei_kernel_return_value; -template struct ei_kernel_impl; +template struct ei_kernel_impl; template struct ei_image_return_value; -template struct ei_image_impl; +template struct ei_image_impl; template class BandMatrix; diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index c5f44dcea..c38789bd9 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -485,21 +485,20 @@ typename ei_traits::Scalar FullPivLU::determinant() cons /********* Implementation of kernel() **************************************************/ -template -struct ei_kernel_impl, Dest> - : ei_kernel_return_value > +template +struct ei_kernel_impl > + : ei_kernel_return_value > { - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; + EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>) + enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime) }; - void evalTo(Dest& dst) const + template void evalTo(Dest& dst) const { - const FullPivLU& dec = this->m_dec; - const int cols = dec.matrixLU().cols(), rank = this->m_rank, dimker = cols - rank; + const int cols = dec().matrixLU().cols(), dimker = cols - rank(); if(dimker == 0) { // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's @@ -525,13 +524,13 @@ struct ei_kernel_impl, Dest> * independent vectors in Ker U. */ - Matrix pivots(rank); - RealScalar premultiplied_threshold = dec.maxPivot() * dec.threshold(); + Matrix pivots(rank()); + RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold(); int p = 0; - for(int i = 0; i < dec.nonzeroPivots(); ++i) - if(ei_abs(dec.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_internal_assert(p == rank); + 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 @@ -539,55 +538,52 @@ struct ei_kernel_impl, Dest> // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified Matrix - m(dec.matrixLU().block(0, 0, rank, cols)); - for(int i = 0; i < rank; ++i) + 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) = dec.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, rank, rank); - m.block(0, 0, rank, rank).template triangularView().setZero(); - for(int i = 0; i < rank; ++i) + m.block(0, 0, rank(), rank()); + m.block(0, 0, rank(), rank()).template triangularView().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, rank, rank) + m.corner(TopLeft, rank(), rank()) .template triangularView().solveInPlace( - m.corner(TopRight, rank, dimker) + m.corner(TopRight, rank(), dimker) ); // now we must undo the column permutation that we had applied! - for(int i = 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 < 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); + 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 -struct ei_image_impl, Dest> - : ei_image_return_value > +template +struct ei_image_impl > + : ei_image_return_value > { - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; + EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>) + enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime) }; - void evalTo(Dest& dst) const + template void evalTo(Dest& dst) const { - const int rank = this->m_rank; - const FullPivLU& dec = this->m_dec; - const MatrixType& originalMatrix = this->m_originalMatrix; - if(rank == 0) + 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 @@ -596,26 +592,28 @@ struct ei_image_impl, Dest> return; } - Matrix pivots(rank); - RealScalar premultiplied_threshold = dec.maxPivot() * dec.threshold(); + Matrix pivots(rank()); + RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold(); int p = 0; - for(int i = 0; i < dec.nonzeroPivots(); ++i) - if(ei_abs(dec.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_internal_assert(p == rank); + ei_internal_assert(p == rank()); - for(int i = 0; i < rank; ++i) - dst.col(i) = originalMatrix.col(dec.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 -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const + EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs) + + template 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: @@ -625,11 +623,9 @@ struct ei_solve_impl, Rhs, Dest> * Step 4: result = Q * c; */ - const FullPivLU& 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 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) @@ -638,35 +634,35 @@ struct ei_solve_impl, Rhs, Dest> return; } - typename Rhs::PlainMatrixType c(rhs.rows(), rhs.cols()); + typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols()); // Step 1 for(int i = 0; i < rows; ++i) - c.row(dec.permutationP().coeff(i)) = rhs.row(i); + c.row(dec().permutationP().coeff(i)) = rhs().row(i); // Step 2 - dec.matrixLU() + dec().matrixLU() .corner(Eigen::TopLeft,smalldim,smalldim) .template triangularView() .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols())); if(rows>cols) { c.corner(Eigen::BottomLeft, rows-cols, c.cols()) - -= dec.matrixLU().corner(Eigen::BottomLeft, rows-cols, cols) + -= dec().matrixLU().corner(Eigen::BottomLeft, rows-cols, cols) * c.corner(Eigen::TopLeft, cols, c.cols()); } // Step 3 - dec.matrixLU() + dec().matrixLU() .corner(TopLeft, nonzero_pivots, nonzero_pivots) .template triangularView() .solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); // Step 4 for(int i = 0; i < nonzero_pivots; ++i) - 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(); + 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/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index eeec3533f..4f6cda68f 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -407,11 +407,13 @@ typename ei_traits::Scalar PartialPivLU::determinant() c /***** Implementation of solve() *****************************************************/ -template -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const + EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs) + + template void evalTo(Dest& dst) const { /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. * So we proceed as follows: @@ -419,23 +421,20 @@ struct ei_solve_impl, Rhs, Dest> * Step 2: replace c by the solution x to Lx = c. * Step 3: replace c by the solution x to Ux = c. */ - - const PartialPivLU& dec = this->m_dec; - const Rhs& rhs = this->m_rhs; - const int size = dec.matrixLU().rows(); - ei_assert(rhs.rows() == size); + const int size = dec().matrixLU().rows(); + ei_assert(rhs().rows() == size); - dst.resize(size, rhs.cols()); + dst.resize(size, rhs().cols()); // Step 1 - for(int i = 0; i < size; ++i) dst.row(dec.permutationP().coeff(i)) = rhs.row(i); + for(int i = 0; i < size; ++i) dst.row(dec().permutationP().coeff(i)) = rhs().row(i); // Step 2 - dec.matrixLU().template triangularView().solveInPlace(dst); + dec().matrixLU().template triangularView().solveInPlace(dst); // Step 3 - dec.matrixLU().template triangularView().solveInPlace(dst); + dec().matrixLU().template triangularView().solveInPlace(dst); } }; diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index a774fdd73..4e2359e0b 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -324,54 +324,52 @@ ColPivHouseholderQR& ColPivHouseholderQR::compute(const return *this; } -template -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const + EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs) + + template void evalTo(Dest& dst) const { - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - const ColPivHouseholderQR& dec = this->m_dec; - const Rhs& rhs = this->m_rhs; - const int rows = dec.rows(), cols = dec.cols(); - dst.resize(cols, rhs.cols()); - ei_assert(rhs.rows() == rows); + const int rows = dec().rows(), cols = dec().cols(); + dst.resize(cols, rhs().cols()); + ei_assert(rhs().rows() == rows); // FIXME introduce nonzeroPivots() and use it here. and more generally, // make the same improvements in this dec as in FullPivLU. - if(dec.rank()==0) + if(dec().rank()==0) { dst.setZero(); return; } - typename Rhs::PlainMatrixType c(rhs); + typename Rhs::PlainMatrixType c(rhs()); // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T c.applyOnTheLeft(makeHouseholderSequence( - dec.matrixQR().corner(TopLeft,rows,dec.rank()), - dec.hCoeffs().start(dec.rank())).transpose() + dec().matrixQR().corner(TopLeft,rows,dec().rank()), + dec().hCoeffs().start(dec().rank())).transpose() ); - if(!dec.isSurjective()) + if(!dec().isSurjective()) { // is c is in the image of R ? - RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff(); - RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec().rank(), c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec().rank(), c.cols()).cwise().abs().maxCoeff(); // FIXME brain dead const RealScalar m_precision = epsilon() * std::min(rows,cols); if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4)) return; } - dec.matrixQR() - .corner(TopLeft, dec.rank(), dec.rank()) + dec().matrixQR() + .corner(TopLeft, dec().rank(), dec().rank()) .template triangularView() - .solveInPlace(c.corner(TopLeft, dec.rank(), c.cols())); + .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); - for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i); - for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero(); + for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i); + for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero(); } }; diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 36ec71b95..ba6866fc0 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -332,57 +332,55 @@ FullPivHouseholderQR& FullPivHouseholderQR::compute(cons return *this; } -template -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const + EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs) + + template void evalTo(Dest& dst) const { - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - const FullPivHouseholderQR& dec = this->m_dec; - const Rhs& rhs = this->m_rhs; - const int rows = dec.rows(), cols = dec.cols(); - dst.resize(cols, rhs.cols()); - ei_assert(rhs.rows() == rows); + const int rows = dec().rows(), cols = dec().cols(); + dst.resize(cols, rhs().cols()); + ei_assert(rhs().rows() == rows); // FIXME introduce nonzeroPivots() and use it here. and more generally, // make the same improvements in this dec as in FullPivLU. - if(dec.rank()==0) + if(dec().rank()==0) { dst.setZero(); return; } - typename Rhs::PlainMatrixType c(rhs); + typename Rhs::PlainMatrixType c(rhs()); - Matrix temp(rhs.cols()); - for (int k = 0; k < dec.rank(); ++k) + Matrix temp(rhs().cols()); + for (int k = 0; k < dec().rank(); ++k) { int remainingSize = rows-k; - c.row(k).swap(c.row(dec.rowsTranspositions().coeff(k))); - c.corner(BottomRight, remainingSize, rhs.cols()) - .applyHouseholderOnTheLeft(dec.matrixQR().col(k).end(remainingSize-1), - dec.hCoeffs().coeff(k), &temp.coeffRef(0)); + c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k))); + c.corner(BottomRight, remainingSize, rhs().cols()) + .applyHouseholderOnTheLeft(dec().matrixQR().col(k).end(remainingSize-1), + dec().hCoeffs().coeff(k), &temp.coeffRef(0)); } - if(!dec.isSurjective()) + if(!dec().isSurjective()) { // is c is in the image of R ? - RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff(); - RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec().rank(), c.cols()).cwise().abs().maxCoeff(); + RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec().rank(), c.cols()).cwise().abs().maxCoeff(); // FIXME brain dead const RealScalar m_precision = epsilon() * std::min(rows,cols); if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision)) return; } - dec.matrixQR() - .corner(TopLeft, dec.rank(), dec.rank()) + dec().matrixQR() + .corner(TopLeft, dec().rank(), dec().rank()) .template triangularView() - .solveInPlace(c.corner(TopLeft, dec.rank(), c.cols())); + .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); - for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i); - for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero(); + for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i); + for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero(); } }; diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 6db0411d9..8742fdf71 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -209,28 +209,28 @@ HouseholderQR& HouseholderQR::compute(const MatrixType& return *this; } -template -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const - { - const HouseholderQR& dec = this->m_dec; - const Rhs& rhs = this->m_rhs; - const int rows = dec.rows(), cols = dec.cols(); - dst.resize(cols, rhs.cols()); - const int rank = std::min(rows, cols); - ei_assert(rhs.rows() == rows); + EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs) - typename Rhs::PlainMatrixType c(rhs); + template void evalTo(Dest& dst) const + { + const int rows = dec().rows(), cols = dec().cols(); + dst.resize(cols, rhs().cols()); + const int rank = std::min(rows, cols); + ei_assert(rhs().rows() == rows); + + typename Rhs::PlainMatrixType c(rhs()); // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T c.applyOnTheLeft(makeHouseholderSequence( - dec.matrixQR().corner(TopLeft,rows,rank), - dec.hCoeffs().start(rank)).transpose() + dec().matrixQR().corner(TopLeft,rows,rank), + dec().hCoeffs().start(rank)).transpose() ); - dec.matrixQR() + dec().matrixQR() .corner(TopLeft, rank, rank) .template triangularView() .solveInPlace(c.corner(TopLeft, rank, c.cols())); diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 8ca425525..2ae7c4859 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -428,33 +428,31 @@ SVD& SVD::compute(const MatrixType& matrix) return *this; } -template -struct ei_solve_impl, Rhs, Dest> - : ei_solve_return_value, Rhs> +template +struct ei_solve_impl, Rhs> + : ei_solve_return_value, Rhs> { - void evalTo(Dest& dst) const + EIGEN_MAKE_SOLVE_HELPERS(SVD<_MatrixType>,Rhs) + + template void evalTo(Dest& dst) const { - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; const int cols = this->cols(); - const SVD& svd = this->m_dec; - const Rhs& rhs = this->m_rhs; - ei_assert(rhs.rows() == svd.rows()); + ei_assert(rhs().rows() == dec().rows()); for (int j=0; j aux = svd.matrixU().adjoint() * rhs.col(j); + Matrix aux = dec().matrixU().adjoint() * rhs().col(j); - for (int i = 0; i struct ei_image_return_value template inline void evalTo(Dest& dst) const { - static_cast *> - (this)->evalTo(dst); + static_cast*>(this)->evalTo(dst); } }; +#define EIGEN_MAKE_IMAGE_HELPERS(DecompositionType) \ + typedef typename DecompositionType::MatrixType MatrixType; \ + typedef typename MatrixType::Scalar Scalar; \ + typedef typename MatrixType::RealScalar RealScalar; \ + inline const DecompositionType& dec() const { return this->m_dec; } \ + inline const MatrixType& originalMatrix() const { return this->m_originalMatrix; } \ + inline int rank() const { return this->m_rank; } + #endif // EIGEN_MISC_IMAGE_H diff --git a/Eigen/src/misc/Kernel.h b/Eigen/src/misc/Kernel.h index bfd75f54b..74ef16abc 100644 --- a/Eigen/src/misc/Kernel.h +++ b/Eigen/src/misc/Kernel.h @@ -63,9 +63,15 @@ template struct ei_kernel_return_value template inline void evalTo(Dest& dst) const { - static_cast *> - (this)->evalTo(dst); + static_cast*>(this)->evalTo(dst); } }; +#define EIGEN_MAKE_KERNEL_HELPERS(DecompositionType) \ + typedef typename DecompositionType::MatrixType MatrixType; \ + typedef typename MatrixType::Scalar Scalar; \ + typedef typename MatrixType::RealScalar RealScalar; \ + inline const DecompositionType& dec() const { return this->m_dec; } \ + inline int rank() const { return this->m_rank; } + #endif // EIGEN_MISC_KERNEL_H diff --git a/Eigen/src/misc/Solve.h b/Eigen/src/misc/Solve.h index c62e34b13..bf8f15773 100644 --- a/Eigen/src/misc/Solve.h +++ b/Eigen/src/misc/Solve.h @@ -57,9 +57,16 @@ template struct ei_solve_return_value template inline void evalTo(Dest& dst) const { - static_cast *> - (this)->evalTo(dst); + static_cast*>(this)->evalTo(dst); } }; +#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType,Rhs) \ + typedef typename DecompositionType::MatrixType MatrixType; \ + typedef typename MatrixType::Scalar Scalar; \ + typedef typename MatrixType::RealScalar RealScalar; \ + typedef typename ei_cleantype::type RhsNestedCleaned; \ + inline const DecompositionType& dec() const { return this->m_dec; } \ + inline const RhsNestedCleaned& rhs() const { return this->m_rhs; } + #endif // EIGEN_MISC_SOLVE_H diff --git a/doc/snippets/FullPivLU_image.cpp b/doc/snippets/FullPivLU_image.cpp index d3092e8b6..817bc1e2d 100644 --- a/doc/snippets/FullPivLU_image.cpp +++ b/doc/snippets/FullPivLU_image.cpp @@ -6,4 +6,4 @@ cout << "Here is the matrix m:" << endl << m << endl; cout << "Notice that the middle column is the sum of the two others, so the " << "columns are linearly dependent." << endl; cout << "Here is a matrix whose columns have the same span but are linearly independent:" - << endl << m.lu().image(m) << endl; + << endl << m.fullPivLu().image(m) << endl; diff --git a/doc/snippets/FullPivLU_kernel.cpp b/doc/snippets/FullPivLU_kernel.cpp index e01186d38..7086e01e2 100644 --- a/doc/snippets/FullPivLU_kernel.cpp +++ b/doc/snippets/FullPivLU_kernel.cpp @@ -1,6 +1,6 @@ MatrixXf m = MatrixXf::Random(3,5); cout << "Here is the matrix m:" << endl << m << endl; -MatrixXf ker = m.lu().kernel(); +MatrixXf ker = m.fullPivLu().kernel(); cout << "Here is a matrix whose columns form a basis of the kernel of m:" << endl << ker << endl; cout << "By definition of the kernel, m*ker is zero:" diff --git a/doc/snippets/FullPivLU_solve.cpp b/doc/snippets/FullPivLU_solve.cpp index ade269789..696f414b3 100644 --- a/doc/snippets/FullPivLU_solve.cpp +++ b/doc/snippets/FullPivLU_solve.cpp @@ -2,7 +2,7 @@ Matrix m = Matrix::Random(); Matrix2f y = Matrix2f::Random(); cout << "Here is the matrix m:" << endl << m << endl; cout << "Here is the matrix y:" << endl << y << endl; -Matrix x = m.lu().solve(y); +Matrix x = m.fillPivLu().solve(y); if((m*x).isApprox(y)) { cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; diff --git a/doc/snippets/HouseholderQR_solve.cpp b/doc/snippets/HouseholderQR_solve.cpp index 429bd81e3..8cce6ce6c 100644 --- a/doc/snippets/HouseholderQR_solve.cpp +++ b/doc/snippets/HouseholderQR_solve.cpp @@ -4,6 +4,6 @@ Matrix3f y = Matrix3f::Random(); cout << "Here is the matrix m:" << endl << m << endl; cout << "Here is the matrix y:" << endl << y << endl; Matrix3f x; -m.householderQr().solve(y, &x); +x = m.householderQr().solve(y); assert(y.isApprox(m*x)); cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;