Closes #1586

Author: dcwhite

src/Core/Algorithms/Math/SolveLinearSystemWithEigen.cc

@@ -42,7 +42,7 @@ using namespace SCIRun::Core;
 
 namespace
 {
-  template <class ColumnMatrixType>
+  template <class ColumnMatrixType, template <typename> class SolverType>
   class SolveLinearSystemAlgorithmEigenCGImpl
   {
   public:
@@ -54,19 +54,19 @@ namespace
     template <class MatrixType>
     typename ColumnMatrixType::EigenBase solveWithEigen(const MatrixType& lhs)
     {
-      Eigen::ConjugateGradient<typename MatrixType::EigenBase> cg;
-      cg.compute(lhs);
+      SolverType<typename MatrixType::EigenBase> solver;
+      solver.compute(lhs);
 
-      if (cg.info() != Eigen::Success)
+      if (solver.info() != Eigen::Success)
         BOOST_THROW_EXCEPTION(AlgorithmInputException()
-          << LinearAlgebraErrorMessage("Conjugate gradient initialization was unsuccessful")
-          << EigenComputationInfo(cg.info()));
-
-      cg.setTolerance(tolerance_);
-      cg.setMaxIterations(maxIterations_);
-      auto solution = cg.solve(*rhs_).eval();
-      tolerance_ = cg.error();
-      maxIterations_ = cg.iterations();
+          << LinearAlgebraErrorMessage("Eigen solver initialization was unsuccessful")
+          << EigenComputationInfo(solver.info()));
+
+      solver.setTolerance(tolerance_);
+      solver.setMaxIterations(maxIterations_);
+      auto solution = solver.solve(*rhs_).eval();
+      tolerance_ = solver.error();
+      maxIterations_ = solver.iterations();
       return solution;
     }
 
@@ -87,6 +87,14 @@ SolveLinearSystemAlgorithm::ComplexOutputs SolveLinearSystemAlgorithm::run(const
   return runImpl<ComplexInputs, ComplexOutputs>(input, params);
 }
 
+template <typename T>
+using CG = Eigen::ConjugateGradient<T>;
+// Not available yet, need to upgrade Eigen
+// template <typename T>
+// using LSCG = Eigen::LeastSquaresConjugateGradient<T>;
+template <typename T>
+using BiCG = Eigen::BiCGSTAB<T>;
+
 template <typename In, typename Out>
 Out SolveLinearSystemAlgorithm::runImpl(const In& input, const Parameters& params) const
 {
@@ -102,8 +110,29 @@ Out SolveLinearSystemAlgorithm::runImpl(const In& input, const Parameters& param
   int maxIterations = std::get<1>(params);
   ENSURE_POSITIVE_INT(maxIterations, "Max iterations out of range!");
 
+  auto method = std::get<2>(params);
+
   using SolutionType = DenseColumnMatrixGeneric<typename std::tuple_element<0, In>::type::element_type::value_type>;
-  using SolverType = SolveLinearSystemAlgorithmEigenCGImpl<SolutionType>;
+  using AlgoTypeCG = SolveLinearSystemAlgorithmEigenCGImpl<SolutionType, CG>;
+  using AlgoTypeBiCG = SolveLinearSystemAlgorithmEigenCGImpl<SolutionType, BiCG>;
+
+  if ("cg" == method)
+    return solve<AlgoTypeCG, In, Out>(input, params);
+  else if ("bicg" == method)
+    return solve<AlgoTypeBiCG, In, Out>(input, params);
+  else
+  {
+    BOOST_THROW_EXCEPTION(AlgorithmProcessingException() << ErrorMessage("Need to upgrade Eigen for LSCG."));
+  }
+}
+
+template <typename SolverType, typename In, typename Out>
+Out SolveLinearSystemAlgorithm::solve(const In& input, const Parameters& params) const
+{
+  auto A = std::get<0>(input);
+  auto b = std::get<1>(input);
+  double tolerance = std::get<0>(params);
+  int maxIterations = std::get<1>(params);
 
   SolverType impl(b, tolerance, maxIterations);
 
@@ -123,8 +152,7 @@ Out SolveLinearSystemAlgorithm::runImpl(const In& input, const Parameters& param
 
   if (x.size() != 0)
   {
-    /// @todo: move ctor
-    auto solution(boost::make_shared<SolutionType>(x));
+    auto solution(boost::make_shared<typename SolverType::SolutionType>(x));
     return Out(solution, impl.tolerance_, impl.maxIterations_);
   }
   else

src/Core/Algorithms/Math/SolveLinearSystemWithEigen.h

@@ -45,7 +45,7 @@ namespace Math {
   public:
     typedef std::tuple<SCIRun::Core::Datatypes::MatrixHandle, SCIRun::Core::Datatypes::DenseColumnMatrixHandle> Inputs;
     typedef std::tuple<SCIRun::Core::Datatypes::ComplexMatrixHandle, SCIRun::Core::Datatypes::ComplexDenseColumnMatrixHandle> ComplexInputs;
-    typedef std::tuple<double, int> Parameters;
+    typedef std::tuple<double, int, std::string> Parameters;
     typedef std::tuple<SCIRun::Core::Datatypes::DenseColumnMatrixHandle, double, int> Outputs;
     typedef std::tuple<SCIRun::Core::Datatypes::ComplexDenseColumnMatrixHandle, double, int> ComplexOutputs;
 
@@ -56,6 +56,8 @@ namespace Math {
   private:
     template <typename In, typename Out>
     Out runImpl(const In& input, const Parameters& params) const;
+    template <typename SolverType, typename In, typename Out>
+    Out solve(const In& input, const Parameters& params) const;
   };
 
   typedef boost::error_info<struct tag_eigen_computation, Eigen::ComputationInfo> EigenComputationInfo;

src/Core/Algorithms/Math/Tests/SolveLinearSystemWithEigenTests.cc

@@ -344,7 +344,7 @@ TEST(EigenSparseSolverTest, DISABLED_CanSolveBigSystem)
     ScopedTimer t("using algorithm object");
     SolveLinearSystemAlgorithm algo;
 
-    x = algo.run(std::make_tuple(A, bCol), std::make_tuple(1e-20, 4000));
+    x = algo.run(std::make_tuple(A, bCol), std::make_tuple(1e-20, 4000, "cg"));
     MatrixHandle solution = std::get<0>(x);
 
     ASSERT_TRUE(solution.get() != nullptr);