[PATCH] solving transposed/adjoint system with SuperLU - now with docs

[Moritz Lenz <moritz@xxxxxxxxxxx>]

---
 Eigen/src/Sparse/SparseLU.h       |   20 +++++++++++++++++---
 Eigen/src/Sparse/SuperLUSupport.h |   18 ++++++++++++++++--
 test/sparse_solvers.cpp           |   10 ++++++++++
 3 files changed, 43 insertions(+), 5 deletions(-)

diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h
index 1425920..46835de 100644
--- a/Eigen/src/Sparse/SparseLU.h
+++ b/Eigen/src/Sparse/SparseLU.h
@@ -25,6 +25,12 @@
 #ifndef EIGEN_SPARSELU_H
 #define EIGEN_SPARSELU_H
 
+enum {
+    SvNoTrans   = 0,
+    SvTranspose = 1,
+    SvAdjoint   = 2
+};
+
 /** \ingroup Sparse_Module
   *
   * \class SparseLU
@@ -115,7 +121,8 @@ class SparseLU
     //inline const MatrixType& matrixU() const { return m_matrixU; }
 
     template<typename BDerived, typename XDerived>
-    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
+    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
+               const int transposed = SvNoTrans) const;
 
     /** \returns true if the factorization succeeded */
     inline bool succeeded(void) const { return m_succeeded; }
@@ -136,10 +143,17 @@ void SparseLU<MatrixType,Backend>::compute(const MatrixType& a)
   ei_assert(false && "not implemented yet");
 }
 
-/** Computes *x = U^-1 L^-1 b */
+/** Computes *x = U^-1 L^-1 b
+ *
+ *  If \transpose transpose is set to SvTranspose or SvAdjoint, the solution
+ *  of the transposed/adjoint system is computed instead.
+ *
+ * Not all backends implement the solution of the transposed or
+ * adjoint system.
+ * */
 template<typename MatrixType, int Backend>
 template<typename BDerived, typename XDerived>
-bool SparseLU<MatrixType,Backend>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const
+bool SparseLU<MatrixType,Backend>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const
 {
   ei_assert(false && "not implemented yet");
   return false;
diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h
index cf0f802..5087a36 100644
--- a/Eigen/src/Sparse/SuperLUSupport.h
+++ b/Eigen/src/Sparse/SuperLUSupport.h
@@ -316,7 +316,7 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
     Scalar determinant() const;
 
     template<typename BDerived, typename XDerived>
-    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
+    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed = SvNoTrans) const;
 
     void compute(const MatrixType& matrix);
 
@@ -411,12 +411,22 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
 
 template<typename MatrixType>
 template<typename BDerived,typename XDerived>
-bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const
+bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b,
+                        MatrixBase<XDerived> *x, const int transposed) const
 {
   const int size = m_matrix.rows();
   const int rhsCols = b.cols();
   ei_assert(size==b.rows());
 
+  switch (transposed) {
+      case SvNoTrans    :  m_sluOptions.Trans = NOTRANS; break;
+      case SvTranspose  :  m_sluOptions.Trans = TRANS;   break;
+      case SvAdjoint    :  m_sluOptions.Trans = CONJ;    break;
+      default:
+        std::cerr << "Eigen: tranpsiotion  option \"" << transposed << "\" not supported by the SuperLU backend\n";
+        m_sluOptions.Trans = NOTRANS;
+  }
+
   m_sluOptions.Fact = FACTORED;
   m_sluOptions.IterRefine = NOREFINE;
 
@@ -441,6 +451,8 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBa
     &m_sluStat, &info, Scalar());
   StatFree(&m_sluStat);
 
+  // reset to previous state
+  m_sluOptions.Trans = NOTRANS;
   return info==0;
 }
 
@@ -563,3 +575,5 @@ typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::dete
 }
 
 #endif // EIGEN_SUPERLUSUPPORT_H
+
+// vim: ts=2 expandtab
diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp
index d1090df..294250e 100644
--- a/test/sparse_solvers.cpp
+++ b/test/sparse_solvers.cpp
@@ -191,6 +191,14 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
           VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU");
         }
         // std::cerr << refDet << " == " << slu.determinant() << "\n";
+        if (slu.solve(b, &x, SvTranspose)) {
+          VERIFY(b.isApprox(m2.transpose() * x, test_precision<Scalar>()));
+        }
+
+        if (slu.solve(b, &x, SvAdjoint)) {
+//          VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>()));
+        }
+
         if (count==0) {
           VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
         }
@@ -227,3 +235,5 @@ void test_sparse_solvers()
     CALL_SUBTEST( sparse_solvers<double>(101, 101) );
   }
 }
+
+// vim: ts=2 sw=2 expandtab
-- 
1.5.6.5


--------------030200080704080904060104--