simplifications in the ei_solve_impl system, factor out some boilerplate code

16 files changed, 187 insertions, 173 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<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
   return *this;
 }
 
-template<typename MatrixType, typename Rhs, typename Dest>
-struct ei_solve_impl<LDLT<MatrixType>, Rhs, Dest>
-  : ei_solve_return_value<LDLT<MatrixType>, Rhs>
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_impl<LDLT<_MatrixType>, Rhs>
+  : ei_solve_return_value<LDLT<_MatrixType>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  EIGEN_MAKE_SOLVE_HELPERS(LDLT<_MatrixType>,Rhs)
+
+  template<typename Dest> 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<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
   return *this;
 }
 
-template<typename MatrixType, int UpLo, typename Rhs, typename Dest>
-struct ei_solve_impl<LLT<MatrixType, UpLo>, Rhs, Dest>
-  : ei_solve_return_value<LLT<MatrixType, UpLo>, Rhs>
+template<typename _MatrixType, int UpLo, typename Rhs>
+struct ei_solve_impl<LLT<_MatrixType, UpLo>, Rhs>
+  : ei_solve_return_value<LLT<_MatrixType, UpLo>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  typedef LLT<_MatrixType,UpLo> LLTType;
+  EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
+
+  template<typename Dest> 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<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, typename Rhs> 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_kernel_impl;
 template<typename DecompositionType> struct ei_image_return_value;
-template<typename DecompositionType, typename Dest> struct ei_image_impl;
+template<typename DecompositionType> 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 c5f44dcea..c38789bd9 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -485,21 +485,20 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons
 
 /********* Implementation of kernel() **************************************************/
 
-template<typename MatrixType, typename Dest>
-struct ei_kernel_impl<FullPivLU<MatrixType>, Dest>
-  : ei_kernel_return_value<FullPivLU<MatrixType> >
+template<typename _MatrixType>
+struct ei_kernel_impl<FullPivLU<_MatrixType> >
+  : ei_kernel_return_value<FullPivLU<_MatrixType> >
 {
-  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<typename Dest> void evalTo(Dest& dst) const
   {
-    const FullPivLU<MatrixType>& 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<FullPivLU<MatrixType>, Dest>
       * independent vectors in Ker U.
       */
 
-    Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank);
-    RealScalar premultiplied_threshold = dec.maxPivot() * dec.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 < 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<FullPivLU<MatrixType>, Dest>
     // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
     Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
            MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
-      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<StrictlyLowerTriangular>().setZero();
-    for(int i = 0; i < 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, rank, rank)
+    m.corner(TopLeft, rank(), rank())
      .template triangularView<UpperTriangular>().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<typename MatrixType, typename Dest>
-struct ei_image_impl<FullPivLU<MatrixType>, Dest>
-  : ei_image_return_value<FullPivLU<MatrixType> >
+template<typename _MatrixType>
+struct ei_image_impl<FullPivLU<_MatrixType> >
+  : ei_image_return_value<FullPivLU<_MatrixType> >
 {
-  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<typename Dest> void evalTo(Dest& dst) const
   {
-    const int rank = this->m_rank;
-    const FullPivLU<MatrixType>& 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<FullPivLU<MatrixType>, Dest>
       return;
     }
     
-    Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank);
-    RealScalar premultiplied_threshold = dec.maxPivot() * dec.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 < 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<typename MatrixType, typename Rhs, typename Dest>
-struct ei_solve_impl<FullPivLU<MatrixType>, Rhs, Dest>
-  : ei_solve_return_value<FullPivLU<MatrixType>, Rhs>
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_impl<FullPivLU<_MatrixType>, Rhs>
+  : ei_solve_return_value<FullPivLU<_MatrixType>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
+  
+  template<typename Dest> 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<FullPivLU<MatrixType>, Rhs, Dest>
      * 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 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<FullPivLU<MatrixType>, 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<UnitLowerTriangular>()
         .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<UpperTriangular>()
         .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<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() c
 
 /***** Implementation of solve() *****************************************************/
 
-template<typename MatrixType, typename Rhs, typename Dest>
-struct ei_solve_impl<PartialPivLU<MatrixType>, Rhs, Dest>
-  : ei_solve_return_value<PartialPivLU<MatrixType>, Rhs>
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_impl<PartialPivLU<_MatrixType>, Rhs>
+  : ei_solve_return_value<PartialPivLU<_MatrixType>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
+  
+  template<typename Dest> 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<PartialPivLU<MatrixType>, 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<MatrixType>& 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<UnitLowerTriangular>().solveInPlace(dst);
+    dec().matrixLU().template triangularView<UnitLowerTriangular>().solveInPlace(dst);
 
     // Step 3
-    dec.matrixLU().template triangularView<UpperTriangular>().solveInPlace(dst);
+    dec().matrixLU().template triangularView<UpperTriangular>().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<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   return *this;
 }
 
-template<typename MatrixType, typename Rhs, typename Dest>
-struct ei_solve_impl<ColPivHouseholderQR<MatrixType>, Rhs, Dest>
-  : ei_solve_return_value<ColPivHouseholderQR<MatrixType>, Rhs>
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_impl<ColPivHouseholderQR<_MatrixType>, Rhs>
+  : ei_solve_return_value<ColPivHouseholderQR<_MatrixType>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
+  
+  template<typename Dest> void evalTo(Dest& dst) const
   {
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    const ColPivHouseholderQR<MatrixType>& 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<Scalar>() * 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<UpperTriangular>()
-       .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<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
   return *this;
 }
 
-template<typename MatrixType, typename Rhs, typename Dest>
-struct ei_solve_impl<FullPivHouseholderQR<MatrixType>, Rhs, Dest>
-  : ei_solve_return_value<FullPivHouseholderQR<MatrixType>, Rhs>
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_impl<FullPivHouseholderQR<_MatrixType>, Rhs>
+  : ei_solve_return_value<FullPivHouseholderQR<_MatrixType>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
   {
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    const FullPivHouseholderQR<MatrixType>& 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<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs.cols());
-    for (int k = 0; k < dec.rank(); ++k)
+    Matrix<Scalar,1,Rhs::ColsAtCompileTime> 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<Scalar>() * 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<UpperTriangular>()
-       .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<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
   return *this;
 }
 
-template<typename MatrixType, typename Rhs, typename Dest>
-struct ei_solve_impl<HouseholderQR<MatrixType>, Rhs, Dest>
-  : ei_solve_return_value<HouseholderQR<MatrixType>, Rhs>
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_impl<HouseholderQR<_MatrixType>, Rhs>
+  : ei_solve_return_value<HouseholderQR<_MatrixType>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
   {
-    const HouseholderQR<MatrixType>& 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 rows = dec().rows(), cols = dec().cols();
+    dst.resize(cols, rhs().cols());
     const int rank = std::min(rows, cols);
-    ei_assert(rhs.rows() == rows);
+    ei_assert(rhs().rows() == rows);
 
-    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,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<UpperTriangular>()
        .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<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
   return *this;
 }
 
-template<typename MatrixType, typename Rhs, typename Dest>
-struct ei_solve_impl<SVD<MatrixType>, Rhs, Dest>
-  : ei_solve_return_value<SVD<MatrixType>, Rhs>
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_impl<SVD<_MatrixType>, Rhs>
+  : ei_solve_return_value<SVD<_MatrixType>, Rhs>
 {
-  void evalTo(Dest& dst) const
+  EIGEN_MAKE_SOLVE_HELPERS(SVD<_MatrixType>,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
   {
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
     const int cols = this->cols();
-    const SVD<MatrixType>& 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<cols; ++j)
     {
-      Matrix<Scalar,MatrixType::RowsAtCompileTime,1> aux = svd.matrixU().adjoint() * rhs.col(j);
+      Matrix<Scalar,MatrixType::RowsAtCompileTime,1> aux = dec().matrixU().adjoint() * rhs().col(j);
 
-      for (int i = 0; i <svd.rows(); ++i)
+      for (int i = 0; i <dec().rows(); ++i)
       {
-        Scalar si = svd.singularValues().coeff(i);
+        Scalar si = dec().singularValues().coeff(i);
         if(si == RealScalar(0))
           aux.coeffRef(i) = Scalar(0);
         else
           aux.coeffRef(i) /= si;
       }
 
-      dst.col(j) = svd.matrixV() * aux;
+      dst.col(j) = dec().matrixV() * aux;
     }
   }
 };
diff --git a/Eigen/src/misc/Image.h b/Eigen/src/misc/Image.h
index a7e2bceec..01e2b160d 100644
--- a/Eigen/src/misc/Image.h
+++ b/Eigen/src/misc/Image.h
@@ -64,9 +64,16 @@ template<typename _DecompositionType> struct ei_image_return_value
 
   template<typename Dest> inline void evalTo(Dest& dst) const
   {
-    static_cast<const ei_image_impl<DecompositionType, Dest> *>
-      (this)->evalTo(dst);
+    static_cast<const ei_image_impl<DecompositionType>*>(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<typename _DecompositionType> struct ei_kernel_return_value
 
   template<typename Dest> inline void evalTo(Dest& dst) const
   {
-    static_cast<const ei_kernel_impl<DecompositionType, Dest> *>
-      (this)->evalTo(dst);
+    static_cast<const ei_kernel_impl<DecompositionType>*>(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<typename _DecompositionType, typename Rhs> struct ei_solve_return_value
 
   template<typename Dest> inline void evalTo(Dest& dst) const
   {
-    static_cast<const ei_solve_impl<DecompositionType, RhsNestedCleaned, Dest> *>
-      (this)->evalTo(dst);
+    static_cast<const ei_solve_impl<DecompositionType,Rhs>*>(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<typename Rhs::Nested>::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<float,2,3> m = Matrix<float,2,3>::Random();
 Matrix2f y = Matrix2f::Random();
 cout << "Here is the matrix m:" << endl << m << endl;
 cout << "Here is the matrix y:" << endl << y << endl;
-Matrix<float,3,2> x = m.lu().solve(y);
+Matrix<float,3,2> 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;