From 8915d5bd2217d4ebbf8fe54ebd2dfbf3f0c420c4 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 1 Feb 2011 17:21:20 +0100
Subject: fix 168 : now TriangularView::solve returns by value making
 TriangularView::solveInPlace less important. Also fix the very outdated
 documentation of this function.

---
 Eigen/src/Core/SolveTriangular.h | 71 ++++++++++++++++++++++++++--------------
 1 file changed, 47 insertions(+), 24 deletions(-)

(limited to 'Eigen/src/Core/SolveTriangular.h')

diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h
index d186b5a4b..b59cdbf85 100644
--- a/Eigen/src/Core/SolveTriangular.h
+++ b/Eigen/src/Core/SolveTriangular.h
@@ -173,8 +173,6 @@ struct triangular_solver_selector<Lhs,Rhs,OnTheRight,Mode,CompleteUnrolling,1> {
 ***************************************************************************/
 
 /** "in-place" version of TriangularView::solve() where the result is written in \a other
-  *
-  *
   *
   * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
   * This function will const_cast it, so constness isn't honored here.
@@ -205,43 +203,68 @@ void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<OtherDerived
 
 /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular.
   *
+  * This function computes the inverse-matrix matrix product inverse(\c *this) * \a other if
+  * \a Side==OnTheLeft (the default), or the right-inverse-multiply  \a other * inverse(\c *this) if
+  * \a Side==OnTheRight.
   *
-  *
-  * This function computes the inverse-matrix matrix product inverse(\c *this) * \a other.
   * The matrix \c *this must be triangular and invertible (i.e., all the coefficients of the
   * diagonal must be non zero). It works as a forward (resp. backward) substitution if \c *this
   * is an upper (resp. lower) triangular matrix.
   *
-  * It is required that \c *this be marked as either an upper or a lower triangular matrix, which
-  * can be done by marked(), and that is automatically the case with expressions such as those returned
-  * by extract().
-  *
   * Example: \include MatrixBase_marked.cpp
   * Output: \verbinclude MatrixBase_marked.out
   *
-  * This function is essentially a wrapper to the faster solveTriangularInPlace() function creating
-  * a temporary copy of \a other, calling solveTriangularInPlace() on the copy and returning it.
-  * Therefore, if \a other is not needed anymore, it is quite faster to call solveTriangularInPlace()
-  * instead of solveTriangular().
+  * This function returns an expression of the inverse-multiply and can works in-place if it is assigned
+  * to the same matrix or vector \a other.
   *
-  * For users coming from BLAS, this function (and more specifically solveTriangularInPlace()) offer
+  * For users coming from BLAS, this function (and more specifically solveInPlace()) offer
   * all the operations supported by the \c *TRSV and \c *TRSM BLAS routines.
   *
-  * \b Tips: to perform a \em "right-inverse-multiply" you can simply transpose the operation, e.g.:
-  * \code
-  * M * T^1  <=>  T.transpose().solveInPlace(M.transpose());
-  * \endcode
-  *
   * \sa TriangularView::solveInPlace()
   */
 template<typename Derived, unsigned int Mode>
-template<int Side, typename RhsDerived>
-typename internal::plain_matrix_type_column_major<RhsDerived>::type
-TriangularView<Derived,Mode>::solve(const MatrixBase<RhsDerived>& rhs) const
+template<int Side, typename Other>
+const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other>
+TriangularView<Derived,Mode>::solve(const MatrixBase<Other>& other) const
 {
-  typename internal::plain_matrix_type_column_major<RhsDerived>::type res(rhs);
-  solveInPlace<Side>(res);
-  return res;
+  return internal::triangular_solve_retval<Side,TriangularView,Other>(*this, other.derived());
 }
 
+namespace internal {
+
+
+template<int Side, typename TriangularType, typename Rhs>
+struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
+{
+  typedef typename internal::plain_matrix_type_column_major<Rhs>::type ReturnType;
+};
+
+template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval
+ : public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
+{
+  typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
+  typedef ReturnByValue<triangular_solve_retval> Base;
+  typedef typename Base::Index Index;
+
+  triangular_solve_retval(const TriangularType& tri, const Rhs& rhs)
+    : m_triangularMatrix(tri), m_rhs(rhs)
+  {}
+
+  inline Index rows() const { return m_rhs.rows(); }
+  inline Index cols() const { return m_rhs.cols(); }
+
+  template<typename Dest> inline void evalTo(Dest& dst) const
+  {
+    if(!(is_same<RhsNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_rhs)))
+      dst = m_rhs;
+    m_triangularMatrix.template solveInPlace<Side>(dst);
+  }
+
+  protected:
+    const TriangularType& m_triangularMatrix;
+    const typename Rhs::Nested m_rhs;
+};
+
+} // namespace internal
+
 #endif // EIGEN_SOLVETRIANGULAR_H
-- 
cgit v1.2.3