trilinos_solver.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#include <deal.II/lac/trilinos_solver.h>

#ifdef DEAL_II_WITH_TRILINOS

#  include <deal.II/base/conditional_ostream.h>
#  include <deal.II/base/std_cxx14/memory.h>

#  include <deal.II/lac/trilinos_precondition.h>
#  include <deal.II/lac/trilinos_sparse_matrix.h>
#  include <deal.II/lac/trilinos_vector.h>

#  include <AztecOO_StatusTest.h>
#  include <AztecOO_StatusTestCombo.h>
#  include <AztecOO_StatusTestMaxIters.h>
#  include <AztecOO_StatusTestResNorm.h>
#  include <AztecOO_StatusType.h>

#  include <cmath>
#  include <limits>

DEAL_II_NAMESPACE_OPEN

namespace TrilinosWrappers
{
  SolverBase::AdditionalData::AdditionalData(
    const bool         output_solver_details,
    const unsigned int gmres_restart_parameter)
    : output_solver_details(output_solver_details)
    , gmres_restart_parameter(gmres_restart_parameter)
  {}



  SolverBase::SolverBase(SolverControl &cn, const AdditionalData &data)
    : solver_name(gmres)
    , solver_control(cn)
    , additional_data(data)
  {}



  SolverBase::SolverBase(const SolverBase::SolverName solver_name,
                         SolverControl &              cn,
                         const AdditionalData &       data)
    : solver_name(solver_name)
    , solver_control(cn)
    , additional_data(data)
  {}



  SolverControl &
  SolverBase::control() const
  {
    return solver_control;
  }



  void
  SolverBase::solve(const SparseMatrix &    A,
                    MPI::Vector &           x,
                    const MPI::Vector &     b,
                    const PreconditionBase &preconditioner)
  {
    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
      &x.trilinos_vector(),
      const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));

    do_solve(preconditioner);
  }



  // Note: "A" is set as a constant reference so that all patterns for ::solve
  //       can be used by the inverse_operator of LinearOperator
  void
  SolverBase::solve(const Epetra_Operator & A,
                    MPI::Vector &           x,
                    const MPI::Vector &     b,
                    const PreconditionBase &preconditioner)
  {
    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_Operator *>(&A),
      &x.trilinos_vector(),
      const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));

    do_solve(preconditioner);
  }



  // Note: "A" is set as a constant reference so that all patterns for ::solve
  //       can be used by the inverse_operator of LinearOperator
  void
  SolverBase::solve(const Epetra_Operator &A,
                    MPI::Vector &          x,
                    const MPI::Vector &    b,
                    const Epetra_Operator &preconditioner)
  {
    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_Operator *>(&A),
      &x.trilinos_vector(),
      const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));

    do_solve(preconditioner);
  }



  // Note: "A" is set as a constant reference so that all patterns for ::solve
  //       can be used by the inverse_operator of LinearOperator
  void
  SolverBase::solve(const Epetra_Operator &   A,
                    Epetra_MultiVector &      x,
                    const Epetra_MultiVector &b,
                    const PreconditionBase &  preconditioner)
  {
    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_Operator *>(&A),
      &x,
      const_cast<Epetra_MultiVector *>(&b));

    do_solve(preconditioner);
  }



  // Note: "A" is set as a constant reference so that all patterns for ::solve
  //       can be used by the inverse_operator of LinearOperator
  void
  SolverBase::solve(const Epetra_Operator &   A,
                    Epetra_MultiVector &      x,
                    const Epetra_MultiVector &b,
                    const Epetra_Operator &   preconditioner)
  {
    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_Operator *>(&A),
      &x,
      const_cast<Epetra_MultiVector *>(&b));

    do_solve(preconditioner);
  }



  void
  SolverBase::solve(const SparseMatrix &          A,
                    ::Vector<double> &      x,
                    const ::Vector<double> &b,
                    const PreconditionBase &      preconditioner)
  {
    // In case we call the solver with deal.II vectors, we create views of the
    // vectors in Epetra format.
    Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
    Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
    Assert(A.local_range().second == A.m(),
           ExcMessage("Can only work in serial when using deal.II vectors."));
    Assert(A.trilinos_matrix().Filled(),
           ExcMessage("Matrix is not compressed. Call compress() method."));

    Epetra_Vector ep_x(View, A.domain_partitioner(), x.begin());
    Epetra_Vector ep_b(View,
                       A.range_partitioner(),
                       const_cast<double *>(b.begin()));

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);

    do_solve(preconditioner);
  }



  void
  SolverBase::solve(Epetra_Operator &             A,
                    ::Vector<double> &      x,
                    const ::Vector<double> &b,
                    const PreconditionBase &      preconditioner)
  {
    Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.begin());
    Epetra_Vector ep_b(View,
                       A.OperatorRangeMap(),
                       const_cast<double *>(b.begin()));

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem =
      std_cxx14::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);

    do_solve(preconditioner);
  }



  void
  SolverBase::solve(const SparseMatrix &                                      A,
                    ::LinearAlgebra::distributed::Vector<double> &      x,
                    const ::LinearAlgebra::distributed::Vector<double> &b,
                    const PreconditionBase &preconditioner)
  {
    // In case we call the solver with deal.II vectors, we create views of the
    // vectors in Epetra format.
    AssertDimension(x.local_size(), A.domain_partitioner().NumMyElements());
    AssertDimension(b.local_size(), A.range_partitioner().NumMyElements());

    Epetra_Vector ep_x(View, A.domain_partitioner(), x.begin());
    Epetra_Vector ep_b(View,
                       A.range_partitioner(),
                       const_cast<double *>(b.begin()));

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);

    do_solve(preconditioner);
  }



  void
  SolverBase::solve(Epetra_Operator &                                         A,
                    ::LinearAlgebra::distributed::Vector<double> &      x,
                    const ::LinearAlgebra::distributed::Vector<double> &b,
                    const PreconditionBase &preconditioner)
  {
    AssertDimension(x.local_size(), A.OperatorDomainMap().NumMyElements());
    AssertDimension(b.local_size(), A.OperatorRangeMap().NumMyElements());

    Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.begin());
    Epetra_Vector ep_b(View,
                       A.OperatorRangeMap(),
                       const_cast<double *>(b.begin()));

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem =
      std_cxx14::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);

    do_solve(preconditioner);
  }


  namespace internal
  {
    namespace
    {
      double
      compute_residual(const Epetra_MultiVector *const residual_vector)
      {
        Assert(residual_vector->NumVectors() == 1,
               ExcMessage("Residual multivector holds more than one vector"));
        TrilinosScalar res_l2_norm = 0.0;
        const int      ierr        = residual_vector->Norm2(&res_l2_norm);
        AssertThrow(ierr == 0, ExcTrilinosError(ierr));
        return res_l2_norm;
      }

      class TrilinosReductionControl : public AztecOO_StatusTest
      {
      public:
        TrilinosReductionControl(const int                   max_steps,
                                 const double                tolerance,
                                 const double                reduction,
                                 const Epetra_LinearProblem &linear_problem);

        virtual ~TrilinosReductionControl() override = default;

        virtual bool
        ResidualVectorRequired() const override
        {
          return status_test_collection->ResidualVectorRequired();
        }

        virtual AztecOO_StatusType
        CheckStatus(int                 CurrentIter,
                    Epetra_MultiVector *CurrentResVector,
                    double              CurrentResNormEst,
                    bool                SolutionUpdated) override
        {
          // Note: CurrentResNormEst is set to -1.0 if no estimate of the
          // residual value is available
          current_residual =
            (CurrentResNormEst < 0.0 ? compute_residual(CurrentResVector) :
                                       CurrentResNormEst);
          if (CurrentIter == 0)
            initial_residual = current_residual;

          return status_test_collection->CheckStatus(CurrentIter,
                                                     CurrentResVector,
                                                     CurrentResNormEst,
                                                     SolutionUpdated);
        }

        virtual AztecOO_StatusType
        GetStatus() const override
        {
          return status_test_collection->GetStatus();
        }

        virtual std::ostream &
        Print(std::ostream &stream, int indent = 0) const override
        {
          return status_test_collection->Print(stream, indent);
        }

        double
        get_initial_residual() const
        {
          return initial_residual;
        }

        double
        get_current_residual() const
        {
          return current_residual;
        }

      private:
        double                                      initial_residual;
        double                                      current_residual;
        std::unique_ptr<AztecOO_StatusTestCombo>    status_test_collection;
        std::unique_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
        std::unique_ptr<AztecOO_StatusTestResNorm>  status_test_abs_tol;
        std::unique_ptr<AztecOO_StatusTestResNorm>  status_test_rel_tol;
      };


      TrilinosReductionControl::TrilinosReductionControl(
        const int                   max_steps,
        const double                tolerance,
        const double                reduction,
        const Epetra_LinearProblem &linear_problem)
        : initial_residual(std::numeric_limits<double>::max())
        , current_residual(std::numeric_limits<double>::max())
        // Consider linear problem converged if any of the collection of
        // criterion are met
        , status_test_collection(
            std_cxx14::make_unique<AztecOO_StatusTestCombo>(
              AztecOO_StatusTestCombo::OR))
      {
        // Maximum number of iterations
        Assert(max_steps >= 0, ExcInternalError());
        status_test_max_steps =
          std_cxx14::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
        status_test_collection->AddStatusTest(*status_test_max_steps);

        Assert(linear_problem.GetRHS()->NumVectors() == 1,
               ExcMessage("RHS multivector holds more than one vector"));

        // Residual norm is below some absolute value
        status_test_abs_tol = std_cxx14::make_unique<AztecOO_StatusTestResNorm>(
          *linear_problem.GetOperator(),
          *(linear_problem.GetLHS()->operator()(0)),
          *(linear_problem.GetRHS()->operator()(0)),
          tolerance);
        status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
                                           AztecOO_StatusTestResNorm::TwoNorm);
        status_test_abs_tol->DefineScaleForm(
          AztecOO_StatusTestResNorm::None, AztecOO_StatusTestResNorm::TwoNorm);
        status_test_collection->AddStatusTest(*status_test_abs_tol);

        // Residual norm, scaled by some initial value, is below some threshold
        status_test_rel_tol = std_cxx14::make_unique<AztecOO_StatusTestResNorm>(
          *linear_problem.GetOperator(),
          *(linear_problem.GetLHS()->operator()(0)),
          *(linear_problem.GetRHS()->operator()(0)),
          reduction);
        status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
                                           AztecOO_StatusTestResNorm::TwoNorm);
        status_test_rel_tol->DefineScaleForm(
          AztecOO_StatusTestResNorm::NormOfInitRes,
          AztecOO_StatusTestResNorm::TwoNorm);
        status_test_collection->AddStatusTest(*status_test_rel_tol);
      }

    } // namespace
  }   // namespace internal


  template <typename Preconditioner>
  void
  SolverBase::do_solve(const Preconditioner &preconditioner)
  {
    int ierr;

    // Next we can allocate the AztecOO solver...
    solver.SetProblem(*linear_problem);

    // ... and we can specify the solver to be used.
    switch (solver_name)
      {
        case cg:
          solver.SetAztecOption(AZ_solver, AZ_cg);
          break;
        case cgs:
          solver.SetAztecOption(AZ_solver, AZ_cgs);
          break;
        case gmres:
          solver.SetAztecOption(AZ_solver, AZ_gmres);
          solver.SetAztecOption(AZ_kspace,
                                additional_data.gmres_restart_parameter);
          break;
        case bicgstab:
          solver.SetAztecOption(AZ_solver, AZ_bicgstab);
          break;
        case tfqmr:
          solver.SetAztecOption(AZ_solver, AZ_tfqmr);
          break;
        default:
          Assert(false, ExcNotImplemented());
      }

    // Set the preconditioner
    set_preconditioner(solver, preconditioner);

    // ... set some options, ...
    solver.SetAztecOption(AZ_output,
                          additional_data.output_solver_details ? AZ_all :
                                                                  AZ_none);
    solver.SetAztecOption(AZ_conv, AZ_noscaled);

    // By default, the Trilinos solver chooses convergence criterion based on
    // the number of iterations made and an absolute tolerance.
    // This implies that the use of the standard Trilinos convergence test
    // actually coincides with ::IterationNumberControl because the
    // solver, unless explicitly told otherwise, will Iterate() until a number
    // of max_steps() are taken or an absolute tolerance() is attained.
    // It is therefore suitable for use with both SolverControl or
    // IterationNumberControl. The final check at the end will determine whether
    // failure to converge to the defined residual norm constitutes failure
    // (SolverControl) or is alright (IterationNumberControl).
    // In the case that the SolverControl wants to perform ReductionControl,
    // then we have to do a little extra something by prescribing a custom
    // status test.
    if (!status_test)
      {
        if (const ReductionControl *const reduction_control =
              dynamic_cast<const ReductionControl *const>(&solver_control))
          {
            status_test =
              std_cxx14::make_unique<internal::TrilinosReductionControl>(
                reduction_control->max_steps(),
                reduction_control->tolerance(),
                reduction_control->reduction(),
                *linear_problem);
            solver.SetStatusTest(status_test.get());
          }
      }

    // ... and then solve!
    ierr =
      solver.Iterate(solver_control.max_steps(), solver_control.tolerance());

    // report errors in more detail than just by checking whether the return
    // status is zero or greater. the error strings are taken from the
    // implementation of the AztecOO::Iterate function
    switch (ierr)
      {
        case -1:
          AssertThrow(false,
                      ExcMessage("AztecOO::Iterate error code -1: "
                                 "option not implemented"));
          break;
        case -2:
          AssertThrow(false,
                      ExcMessage("AztecOO::Iterate error code -2: "
                                 "numerical breakdown"));
          break;
        case -3:
          AssertThrow(false,
                      ExcMessage("AztecOO::Iterate error code -3: "
                                 "loss of precision"));
          break;
        case -4:
          AssertThrow(false,
                      ExcMessage("AztecOO::Iterate error code -4: "
                                 "GMRES Hessenberg ill-conditioned"));
          break;
        default:
          AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
      }

    // Finally, let the deal.II SolverControl object know what has
    // happened. If the solve succeeded, the status of the solver control will
    // turn into SolverControl::success.
    // If the residual is not computed/stored by the solver, as can happen for
    // certain choices of solver or if a custom status test is set, then the
    // result returned by TrueResidual() is equal to -1. In this case we must
    // compute it ourself.
    if (const internal::TrilinosReductionControl
          *const reduction_control_status =
            dynamic_cast<const internal::TrilinosReductionControl *const>(
              status_test.get()))
      {
        Assert(dynamic_cast<const ReductionControl *const>(&solver_control),
               ExcInternalError());

        // Check to see if solver converged in one step
        // This can happen if the matrix is diagonal and a non-trivial
        // preconditioner is used.
        if (solver.NumIters() > 0)
          {
            // For ReductionControl, we must first register the initial residual
            // value. This is the basis from which it will determine whether the
            // current residual corresponds to a converged state.
            solver_control.check(
              0, reduction_control_status->get_initial_residual());
            solver_control.check(
              solver.NumIters(),
              reduction_control_status->get_current_residual());
          }
        else
          solver_control.check(
            solver.NumIters(),
            reduction_control_status->get_current_residual());
      }
    else
      {
        Assert(solver.TrueResidual() >= 0.0, ExcInternalError());
        solver_control.check(solver.NumIters(), solver.TrueResidual());
      }

    if (solver_control.last_check() != SolverControl::success)
      AssertThrow(false,
                  SolverControl::NoConvergence(solver_control.last_step(),
                                               solver_control.last_value()));
  }



  template <>
  void
  SolverBase::set_preconditioner(AztecOO &               solver,
                                 const PreconditionBase &preconditioner)
  {
    // Introduce the preconditioner, if the identity preconditioner is used,
    // the precondioner is set to none, ...
    if (preconditioner.preconditioner.strong_count() != 0)
      {
        const int ierr = solver.SetPrecOperator(
          const_cast<Epetra_Operator *>(preconditioner.preconditioner.get()));
        AssertThrow(ierr == 0, ExcTrilinosError(ierr));
      }
    else
      solver.SetAztecOption(AZ_precond, AZ_none);
  }


  template <>
  void
  SolverBase::set_preconditioner(AztecOO &              solver,
                                 const Epetra_Operator &preconditioner)
  {
    const int ierr =
      solver.SetPrecOperator(const_cast<Epetra_Operator *>(&preconditioner));
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
  }


  /* ---------------------- SolverCG ------------------------ */

  SolverCG::AdditionalData::AdditionalData(const bool output_solver_details)
    : SolverBase::AdditionalData(output_solver_details)
  {}



  SolverCG::SolverCG(SolverControl &cn, const AdditionalData &data)
    : SolverBase(cn, data)
    , additional_data(data)
  {
    solver_name = cg;
  }


  /* ---------------------- SolverGMRES ------------------------ */

  SolverGMRES::AdditionalData::AdditionalData(
    const bool         output_solver_details,
    const unsigned int restart_parameter)
    : SolverBase::AdditionalData(output_solver_details, restart_parameter)
  {}



  SolverGMRES::SolverGMRES(SolverControl &cn, const AdditionalData &data)
    : SolverBase(cn, data)
    , additional_data(data)
  {
    solver_name = gmres;
  }


  /* ---------------------- SolverBicgstab ------------------------ */

  SolverBicgstab::AdditionalData::AdditionalData(
    const bool output_solver_details)
    : SolverBase::AdditionalData(output_solver_details)
  {}



  SolverBicgstab::SolverBicgstab(SolverControl &cn, const AdditionalData &data)
    : SolverBase(cn, data)
    , additional_data(data)
  {
    solver_name = bicgstab;
  }


  /* ---------------------- SolverCGS ------------------------ */

  SolverCGS::AdditionalData::AdditionalData(const bool output_solver_details)
    : SolverBase::AdditionalData(output_solver_details)
  {}



  SolverCGS::SolverCGS(SolverControl &cn, const AdditionalData &data)
    : SolverBase(cn, data)
    , additional_data(data)
  {
    solver_name = cgs;
  }


  /* ---------------------- SolverTFQMR ------------------------ */

  SolverTFQMR::AdditionalData::AdditionalData(const bool output_solver_details)
    : SolverBase::AdditionalData(output_solver_details)
  {}



  SolverTFQMR::SolverTFQMR(SolverControl &cn, const AdditionalData &data)
    : SolverBase(cn, data)
    , additional_data(data)
  {
    solver_name = tfqmr;
  }



  /* ---------------------- SolverDirect ------------------------ */

  SolverDirect::AdditionalData::AdditionalData(const bool output_solver_details,
                                               const std::string &solver_type)
    : output_solver_details(output_solver_details)
    , solver_type(solver_type)
  {}



  SolverDirect::SolverDirect(SolverControl &cn, const AdditionalData &data)
    : solver_control(cn)
    , additional_data(data.output_solver_details, data.solver_type)
  {}



  SolverControl &
  SolverDirect::control() const
  {
    return solver_control;
  }



  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 = std_cxx14::make_unique<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;

    // First set whether we want to print the solver information to screen or
    // not.
    ConditionalOStream verbose_cout(std::cout,
                                    additional_data.output_solver_details);

    // 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));
  }


  void
  SolverDirect::solve(MPI::Vector &x, const MPI::Vector &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()));

    // 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 = solver->Solve();
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));

    // Finally, force the SolverControl object to report convergence
    solver_control.check(0, 0);
  }



  void
  SolverDirect::solve(
    ::LinearAlgebra::distributed::Vector<double> &      x,
    const ::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;

    // Fetch return value of Amesos Solver functions
    int 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);

    // 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();
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));

    // Finally, let the deal.II SolverControl object know what has
    // happened. If the solve succeeded, the status of the solver control will
    // turn into SolverControl::success.
    solver_control.check(0, 0);

    if (solver_control.last_check() != SolverControl::success)
      AssertThrow(false,
                  SolverControl::NoConvergence(solver_control.last_step(),
                                               solver_control.last_value()));
  }


  void
  SolverDirect::solve(const SparseMatrix &A,
                      MPI::Vector &       x,
                      const MPI::Vector & b)
  {
    // We need an Epetra_LinearProblem object to let the Amesos solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
      &x.trilinos_vector(),
      const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));

    do_solve();
  }



  void
  SolverDirect::solve(const SparseMatrix &          A,
                      ::Vector<double> &      x,
                      const ::Vector<double> &b)
  {
    // In case we call the solver with deal.II vectors, we create views of the
    // vectors in Epetra format.
    Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
    Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
    Assert(A.local_range().second == A.m(),
           ExcMessage("Can only work in serial when using deal.II vectors."));
    Epetra_Vector ep_x(View, A.domain_partitioner(), x.begin());
    Epetra_Vector ep_b(View,
                       A.range_partitioner(),
                       const_cast<double *>(b.begin()));

    // We need an Epetra_LinearProblem object to let the Amesos solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);

    do_solve();
  }



  void
  SolverDirect::solve(
    const SparseMatrix &                                      A,
    ::LinearAlgebra::distributed::Vector<double> &      x,
    const ::LinearAlgebra::distributed::Vector<double> &b)
  {
    AssertDimension(x.local_size(), A.domain_partitioner().NumMyElements());
    AssertDimension(b.local_size(), A.range_partitioner().NumMyElements());
    Epetra_Vector ep_x(View, A.domain_partitioner(), x.begin());
    Epetra_Vector ep_b(View,
                       A.range_partitioner(),
                       const_cast<double *>(b.begin()));

    // We need an Epetra_LinearProblem object to let the Amesos solver know
    // about the matrix and vectors.
    linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
      const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);

    do_solve();
  }
} // namespace TrilinosWrappers


// explicit instantiations
// TODO: put these instantiations into generic file
namespace TrilinosWrappers
{
  template void
  SolverBase::do_solve(const PreconditionBase &preconditioner);

  template void
  SolverBase::do_solve(const Epetra_Operator &preconditioner);
} // namespace TrilinosWrappers

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_PETSC