Add internal method _solve_impl_transposed() to LU decomposition classes that solves A^T x = b or A^* x = b.

1 files changed, 33 insertions, 3 deletions

diff --git a/test/lu.cpp b/test/lu.cpp
index b90367438..f753fc74a 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -92,6 +92,20 @@ template<typename MatrixType> void lu_non_invertible()
   // test that the code, which does resize(), may be applied to an xpr
   m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3);
   VERIFY_IS_APPROX(m3, m1*m2);
+
+  // test solve with transposed
+  m3 = MatrixType::Random(rows,cols2);
+  m2 = m1.transpose()*m3;
+  m3 = MatrixType::Random(rows,cols2);
+  lu.template _solve_impl_transposed<false>(m2, m3);
+  VERIFY_IS_APPROX(m2, m1.transpose()*m3);
+
+  // test solve with conjugate transposed
+  m3 = MatrixType::Random(rows,cols2);
+  m2 = m1.adjoint()*m3;
+  m3 = MatrixType::Random(rows,cols2);
+  lu.template _solve_impl_transposed<true>(m2, m3);
+  VERIFY_IS_APPROX(m2, m1.adjoint()*m3);
 }
 
 template<typename MatrixType> void lu_invertible()
@@ -124,6 +138,12 @@ template<typename MatrixType> void lu_invertible()
   m2 = lu.solve(m3);
   VERIFY_IS_APPROX(m3, m1*m2);
   VERIFY_IS_APPROX(m2, lu.inverse()*m3);
+  // test solve with transposed
+  lu.template _solve_impl_transposed<false>(m3, m2);
+  VERIFY_IS_APPROX(m3, m1.transpose()*m2);
+  // test solve with conjugate transposed
+  lu.template _solve_impl_transposed<true>(m3, m2);
+  VERIFY_IS_APPROX(m3, m1.adjoint()*m2);
 
   // Regression test for Bug 302
   MatrixType m4 = MatrixType::Random(size,size);
@@ -136,14 +156,24 @@ template<typename MatrixType> void lu_partial_piv()
      PartialPivLU.h
   */
   typedef typename MatrixType::Index Index;
-  Index rows = internal::random<Index>(1,4);
-  Index cols = rows;
+  Index size = internal::random<Index>(1,4);
 
-  MatrixType m1(cols, rows);
+  MatrixType m1(size, size), m2(size, size), m3(size, size);
   m1.setRandom();
   PartialPivLU<MatrixType> plu(m1);
 
   VERIFY_IS_APPROX(m1, plu.reconstructedMatrix());
+
+  m3 = MatrixType::Random(size,size);
+  m2 = plu.solve(m3);
+  VERIFY_IS_APPROX(m3, m1*m2);
+  VERIFY_IS_APPROX(m2, plu.inverse()*m3);
+  // test solve with transposed
+  plu.template _solve_impl_transposed<false>(m3, m2);
+  VERIFY_IS_APPROX(m3, m1.transpose()*m2);
+  // test solve with conjugate transposed
+  plu.template _solve_impl_transposed<true>(m3, m2);
+  VERIFY_IS_APPROX(m3, m1.adjoint()*m2);
 }
 
 template<typename MatrixType> void lu_verify_assert()