From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Fri, 10 Nov 2017 13:50:09 +0000 (+0100)
Subject: Enable Trilinos direct solvers with LA::distributed::Vector.
X-Git-Tag: v9.0.0-rc1~790^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=784caee7cb8d25ed44967140ff6cafec9bca5199;p=dealii.git

Enable Trilinos direct solvers with LA::distributed::Vector.
---

diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h
index 17e672d342..38d63ac5de 100644
--- a/include/deal.II/lac/trilinos_solver.h
+++ b/include/deal.II/lac/trilinos_solver.h
@@ -645,6 +645,15 @@ namespace TrilinosWrappers
      */
     void solve (MPI::Vector &x, const MPI::Vector &b);
 
+    /**
+     * Solve the linear system <tt>Ax=b</tt> based on the package set in
+     * intialize() for deal.II's own parallel vectors. Note the matrix is not
+     * refactorized during this call.
+     */
+    void
+    solve (dealii::LinearAlgebra::distributed::Vector<double>       &x,
+           const dealii::LinearAlgebra::distributed::Vector<double> &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 2839e10dc2..79c19a658e 100644
--- a/source/lac/trilinos_solver.cc
+++ b/source/lac/trilinos_solver.cc
@@ -780,17 +780,45 @@ namespace TrilinosWrappers
     // Assign the RHS vector to the Epetra_LinearProblem object
     linear_problem->SetRHS(const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
 
+    // 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;
+
     // Fetch return value of Amesos Solver functions
-    int ierr;
+    int ierr = solver->Solve ();
+    AssertThrow (ierr == 0, ExcTrilinosError (ierr));
+
+    // Finally, force the SolverControl object to report convergence
+    solver_control.check (0, 0);
+  }
+
+
+
+  void SolverDirect::solve (dealii::LinearAlgebra::distributed::Vector<double>       &x,
+                            const dealii::LinearAlgebra::distributed::Vector<double> &b)
+  {
+    Epetra_Vector ep_x (View, linear_problem->GetOperator()->OperatorDomainMap(), x.begin());
+    Epetra_Vector ep_b (View, linear_problem->GetOperator()->OperatorRangeMap(), const_cast<double *>(b.begin()));
+
+    // Assign the empty LHS vector to the Epetra_LinearProblem object
+    linear_problem->SetLHS(&ep_x);
+
+    // Assign the RHS vector to the Epetra_LinearProblem object
+    linear_problem->SetRHS(&ep_b);
 
     // 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 ();
+
+    // Fetch return value of Amesos Solver functions
+    int ierr = solver->Solve ();
     AssertThrow (ierr == 0, ExcTrilinosError (ierr));
 
     // Finally, force the SolverControl object to report convergence
@@ -917,7 +945,6 @@ namespace TrilinosWrappers
 
     do_solve();
   }
-
 }