From 0e9866f59fe84127a8ca5216a0e3255d747e41ed Mon Sep 17 00:00:00 2001
From: Mike Harmon <mdh266@gmail.com>
Date: Sun, 10 Apr 2016 16:54:19 +0200
Subject: [PATCH] Refactored trilinos direct solver

---
 include/deal.II/lac/trilinos_solver.h | 15 ++++++
 source/lac/trilinos_solver.cc         | 77 ++++++++++++++++++++++++++-
 2 files changed, 90 insertions(+), 2 deletions(-)

diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h
index 9ab13fc2aa..084bf77a41 100644
--- a/include/deal.II/lac/trilinos_solver.h
+++ b/include/deal.II/lac/trilinos_solver.h
@@ -567,6 +567,21 @@ namespace TrilinosWrappers
      */
     virtual ~SolverDirect ();
 
+    /**
+     * Initializes the direct solver for the matrix <tt>A</tt> and creates a
+     * factorization for it with the package chosen from the additional
+     * data structure. Note that there is no need for a preconditioner
+     * here and solve() is not called.
+     */
+    void initialize (const SparseMatrix &A);
+
+    /*
+     * Solve the linear system <tt>Ax=b</tt> based on the
+     * package set in intialize(). Note the matrix is not refactorized during
+     * this call.
+     */
+    void solve (VectorBase &x, const VectorBase &b);
+
     /**
      * Solve the linear system <tt>Ax=b</tt>. Creates a factorization of the
      * matrix with the package chosen from the additional data structure and
diff --git a/source/lac/trilinos_solver.cc b/source/lac/trilinos_solver.cc
index 6b539705ce..b28637bc55 100644
--- a/source/lac/trilinos_solver.cc
+++ b/source/lac/trilinos_solver.cc
@@ -442,9 +442,15 @@ namespace TrilinosWrappers
 
 
 
-  void
-  SolverDirect::do_solve()
+  void SolverDirect::initialize (const SparseMatrix &A)
   {
+    // We need an Epetra_LinearProblem object to let the Amesos solver know
+    // about the matrix and vectors.
+    linear_problem.reset (new Epetra_LinearProblem ());
+
+    // Assign the matrix operator to the Epetra_LinearProblem object
+    linear_problem->SetOperator(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()));
+
     // Fetch return value of Amesos Solver functions
     int ierr;
 
@@ -480,6 +486,73 @@ namespace TrilinosWrappers
     verbose_cout << "Starting numeric factorization" << std::endl;
     ierr = solver->NumericFactorization();
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+  void SolverDirect::solve (VectorBase &x, const VectorBase &b)
+  {
+    // Assign the empty LHS vector to the Epetra_LinearProblem object
+    linear_problem->SetLHS(&x.trilinos_vector());
+
+    // Assign the RHS vector to the Epetra_LinearProblem object
+    linear_problem->SetRHS(const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
+
+    // Fetch return value of Amesos Solver functions
+    int ierr;
+
+    // First set whether we want to print the solver information to screen or
+    // not.
+    ConditionalOStream verbose_cout (std::cout,
+                                     additional_data.output_solver_details);
+
+
+    verbose_cout << "Starting solve" << std::endl;
+    ierr = solver->Solve ();
+    AssertThrow (ierr == 0, ExcTrilinosError (ierr));
+
+    // Finally, force the SolverControl object to report convergence
+    solver_control.check (0, 0);
+  }
+
+
+
+  void
+  SolverDirect::do_solve()
+  {
+    // Fetch return value of Amesos Solver functions
+    int ierr;
+
+    // First set whether we want to print the solver information to screen or
+    // not.
+    ConditionalOStream verbose_cout (std::cout,
+                                     additional_data.output_solver_details);
+
+    solver.reset ();
+
+    // Next allocate the Amesos solver, this is done in two steps, first we
+    // create a solver Factory and and generate with that the concrete Amesos
+    // solver, if possible.
+    Amesos Factory;
+
+    AssertThrow (Factory.Query (additional_data.solver_type.c_str ()),
+                 ExcMessage (std::
+                             string ("You tried to select the solver type <")
+                             + additional_data.solver_type +
+                             "> but this solver is not supported by Trilinos either "
+                             "because it does not exist, or because Trilinos was not "
+                             "configured for its use."));
+
+    solver.reset (Factory.
+                  Create (additional_data.solver_type.c_str (),
+                          *linear_problem));
+
+    verbose_cout << "Starting symbolic factorization" << std::endl;
+    ierr = solver->SymbolicFactorization ();
+    AssertThrow (ierr == 0, ExcTrilinosError (ierr));
+
+    verbose_cout << "Starting numeric factorization" << std::endl;
+    ierr = solver->NumericFactorization ();
+    AssertThrow (ierr == 0, ExcTrilinosError (ierr));
 
     verbose_cout << "Starting solve" << std::endl;
     ierr = solver->Solve();
-- 
2.39.5