trilinos_solver.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

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

#ifdef DEAL_II_WITH_TRILINOS

#  include <deal.II/base/conditional_ostream.h>
#  include <deal.II/lac/trilinos_sparse_matrix.h>
#  include <deal.II/lac/trilinos_vector_base.h>
#  include <deal.II/lac/trilinos_precondition.h>

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

#  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)
    :
    solver_name    (gmres),
    solver_control (cn)
  {}



  SolverBase::SolverBase (const enum SolverBase::SolverName  solver_name,
                          SolverControl                     &cn)
    :
    solver_name    (solver_name),
    solver_control (cn)
  {}



  SolverBase::~SolverBase ()
  {}



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



  void
  SolverBase::solve (const SparseMatrix     &A,
                     VectorBase             &x,
                     const VectorBase       &b,
                     const PreconditionBase &preconditioner)
  {
    linear_problem.reset();

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem.reset
    (new 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,
                     VectorBase             &x,
                     const VectorBase       &b,
                     const PreconditionBase &preconditioner)
  {
    linear_problem.reset();

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem.reset
    (new 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,
                     VectorBase             &x,
                     const VectorBase       &b,
                     const Epetra_Operator  &preconditioner)
  {
    linear_problem.reset();

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem.reset
    (new 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)
  {
    linear_problem.reset();

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem.reset
    (new 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)
  {
    linear_problem.reset();

    // We need an Epetra_LinearProblem object to let the AztecOO solver know
    // about the matrix and vectors.
    linear_problem.reset
    (new 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)
  {
    linear_problem.reset();

    // 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.reset (new 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)
  {
    linear_problem.reset();

    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.reset (new 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)
  {
    linear_problem.reset();

    // In case we call the solver with deal.II vectors, we create views of the
    // vectors in Epetra format.
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
                     A.domain_partitioner().NumMyElements());
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(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.reset (new 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)
  {
    linear_problem.reset();

    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
                     A.OperatorDomainMap().NumMyElements());
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(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.reset (new 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 double               &max_steps,
                                  const double               &tolerance,
                                  const double               &reduction,
                                  const Epetra_LinearProblem &linear_problem);

        virtual ~TrilinosReductionControl() {}

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

        virtual AztecOO_StatusType
        CheckStatus (int                 CurrentIter,
                     Epetra_MultiVector *CurrentResVector,
                     double              CurrentResNormEst,
                     bool                SolutionUpdated)
        {
          // 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
        {
          return status_test_collection->GetStatus();
        }

        virtual std::ostream &
        Print (std::ostream &stream,
               int           indent = 0) const
        {
          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_cxx11::shared_ptr<AztecOO_StatusTestCombo>    status_test_collection;
        std_cxx11::shared_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
        std_cxx11::shared_ptr<AztecOO_StatusTestResNorm>  status_test_abs_tol;
        std_cxx11::shared_ptr<AztecOO_StatusTestResNorm>  status_test_rel_tol;
      };


      TrilinosReductionControl::TrilinosReductionControl(
        const double               &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.reset(
          new AztecOO_StatusTestCombo (AztecOO_StatusTestCombo::OR) );

        // Maximum number of iterations
        status_test_max_steps.reset(
          new 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.reset(
          new 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.reset(
          new 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);
      }

    }
  }


  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.reset(new 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"));
      case -2:
        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -2: "
                                       "numerical breakdown"));
      case -3:
        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -3: "
                                       "loss of precision"));
      case -4:
        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -4: "
                                       "GMRES Hessenberg ill-conditioned"));
      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.use_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)
    :
    output_solver_details (output_solver_details)
  {}



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


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

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



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


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

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




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


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

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




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


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

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



  SolverTFQMR::SolverTFQMR (SolverControl        &cn,
                            const AdditionalData &data)
    :
    SolverBase (cn),
    additional_data (data.output_solver_details)
  {
    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)
  {}



  SolverDirect::~SolverDirect ()
  {}



  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.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;

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


  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();
    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,
                       VectorBase             &x,
                       const VectorBase       &b)
  {
    // We need an Epetra_LinearProblem object to let the Amesos solver know
    // about the matrix and vectors.
    linear_problem.reset
    (new 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.reset (new 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 (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
                     A.domain_partitioner().NumMyElements());
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(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.reset (new Epetra_LinearProblem
                          (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
                           &ep_x, &ep_b));

    do_solve();
  }

}


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

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_PETSC