patch from Moritz Lenz to allow solving transposed problem with superlu

3 files changed, 39 insertions, 5 deletions

diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h
index 142592050..9cab67a51 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 \a 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 6df94e35b..2aeccc534 100644
--- a/Eigen/src/Sparse/SuperLUSupport.h
+++ b/Eigen/src/Sparse/SuperLUSupport.h
@@ -318,7 +318,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);
 
@@ -413,12 +413,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;
 
@@ -443,6 +453,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;
 }
 
diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp
index d1090dfed..e1ec1ef35 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
         }