[PATCH] transposed/conjugated solving for SuperLU
• Subject: [PATCH] transposed/conjugated solving for SuperLU
• From: Moritz Lenz <moritz@xxxxxxxxxxx>
• Date: Mon, 6 Apr 2009 10:55:47 +0200

```---
Eigen/src/Sparse/SparseLU.h       |   11 +++++++++--
Eigen/src/Sparse/SuperLUSupport.h |   18 ++++++++++++++++--
test/sparse_solvers.cpp           |   10 ++++++++++
3 files changed, 35 insertions(+), 4 deletions(-)

diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h
index 1425920..80f8490 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,
+};
+
/** \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; }
@@ -139,7 +146,7 @@ void SparseLU<MatrixType,Backend>::compute(const MatrixType& a)
/** Computes *x = U^-1 L^-1 b */
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)) {
+        }
+
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

--------------040703010901090908020503--

```

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/