trilinos_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2016 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.
//
// ---------------------------------------------------------------------

#ifndef dealii__trilinos_vector_h
#define dealii__trilinos_vector_h


#include <deal.II/base/config.h>

#ifdef DEAL_II_WITH_TRILINOS

#  include <deal.II/base/std_cxx11/shared_ptr.h>
#  include <deal.II/base/subscriptor.h>
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/utilities.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/trilinos_vector_base.h>

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#  include "Epetra_Map.h"
#  include "Epetra_LocalMap.h"
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

DEAL_II_NAMESPACE_OPEN


// forward declaration
template <typename> class Vector;

namespace TrilinosWrappers
{
  class SparseMatrix;

  namespace
  {
#ifndef DEAL_II_WITH_64BIT_INDICES
    // define a helper function that queries the global ID of local ID of
    // an Epetra_BlockMap object  by calling either the 32- or 64-bit
    // function necessary.
    inline
    int gid(const Epetra_BlockMap &map, int i)
    {
      return map.GID(i);
    }
#else
    // define a helper function that queries the global ID of local ID of
    // an Epetra_BlockMap object  by calling either the 32- or 64-bit
    // function necessary.
    inline
    long long int gid(const Epetra_BlockMap &map, int i)
    {
      return map.GID64(i);
    }
#endif
  }

  namespace MPI
  {
    class BlockVector;

    class Vector : public VectorBase
    {
    public:
      typedef ::types::global_dof_index size_type;

      static const bool supports_distributed_data = true;

      Vector ();

      Vector (const Vector &v);

#ifdef DEAL_II_WITH_CXX11

      Vector (Vector &&v);
#endif

      ~Vector ();

      void reinit (const VectorBase &v,
                   const bool        omit_zeroing_entries = false,
                   const bool        allow_different_maps = false);

      void reinit (const BlockVector &v,
                   const bool         import_data = false);

      Vector &operator= (const TrilinosScalar s);

      Vector &operator= (const Vector &v);

#ifdef DEAL_II_WITH_CXX11

      Vector &operator= (Vector &&v);
#endif

      Vector &operator= (const ::TrilinosWrappers::Vector &v);

      template <typename Number>
      Vector &operator= (const ::Vector<Number> &v);

      void import_nonlocal_data_for_fe
      (const ::TrilinosWrappers::SparseMatrix &matrix,
       const Vector                                 &vector);

      explicit Vector (const Epetra_Map &parallel_partitioning)  DEAL_II_DEPRECATED;

      Vector (const Epetra_Map &parallel_partitioning,
              const VectorBase &v) DEAL_II_DEPRECATED;

      template <typename number>
      void reinit (const Epetra_Map             &parallel_partitioner,
                   const ::Vector<number> &v) DEAL_II_DEPRECATED;

      void reinit (const Epetra_Map &parallel_partitioning,
                   const bool        omit_zeroing_entries = false) DEAL_II_DEPRECATED;

      template <typename Number>
      Vector (const Epetra_Map             &parallel_partitioning,
              const ::Vector<Number> &v) DEAL_II_DEPRECATED;

      explicit Vector (const IndexSet &parallel_partitioning,
                       const MPI_Comm &communicator = MPI_COMM_WORLD);

      Vector (const IndexSet &local,
              const IndexSet &ghost,
              const MPI_Comm &communicator = MPI_COMM_WORLD);

      Vector (const IndexSet   &parallel_partitioning,
              const VectorBase &v,
              const MPI_Comm   &communicator = MPI_COMM_WORLD);

      template <typename Number>
      Vector (const IndexSet               &parallel_partitioning,
              const ::Vector<Number> &v,
              const MPI_Comm               &communicator = MPI_COMM_WORLD);

      void reinit (const IndexSet &parallel_partitioning,
                   const MPI_Comm &communicator = MPI_COMM_WORLD,
                   const bool      omit_zeroing_entries = false);

      void reinit (const IndexSet &locally_owned_entries,
                   const IndexSet &ghost_entries,
                   const MPI_Comm &communicator = MPI_COMM_WORLD,
                   const bool      vector_writable = false);
    };




// ------------------- inline and template functions --------------


    inline
    void swap (Vector &u, Vector &v)
    {
      u.swap (v);
    }


#ifndef DOXYGEN

    template <typename number>
    Vector::Vector (const Epetra_Map             &input_map,
                    const ::Vector<number> &v)
    {
      reinit (input_map, v);
    }



    template <typename number>
    Vector::Vector (const IndexSet               &parallel_partitioner,
                    const ::Vector<number> &v,
                    const MPI_Comm               &communicator)
    {
      *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
                     v);
    }




    template <typename number>
    void Vector::reinit (const Epetra_Map             &parallel_partitioner,
                         const ::Vector<number> &v)
    {
      if (vector.get() == 0 || vector->Map().SameAs(parallel_partitioner) == false)
        vector.reset (new Epetra_FEVector(parallel_partitioner));

      has_ghosts = vector->Map().UniqueGIDs()==false;

      const int size = parallel_partitioner.NumMyElements();

      // Need to copy out values, since the deal.II might not use doubles, so
      // that a direct access is not possible.
      for (int i=0; i<size; ++i)
        (*vector)[0][i] = v(gid(parallel_partitioner,i));
    }


    inline
    Vector &
    Vector::operator= (const TrilinosScalar s)
    {
      VectorBase::operator= (s);

      return *this;
    }


    template <typename Number>
    Vector &
    Vector::operator= (const ::Vector<Number> &v)
    {
      if (size() != v.size())
        {
          vector.reset (new Epetra_FEVector(Epetra_Map
                                            (static_cast<TrilinosWrappers::types::int_type>(v.size()), 0,
#ifdef DEAL_II_WITH_MPI
                                             Epetra_MpiComm(MPI_COMM_SELF)
#else
                                             Epetra_SerialComm()
#endif
                                            )));
        }

      reinit (vector_partitioner(), v);
      return *this;
    }


#endif /* DOXYGEN */

  } /* end of namespace MPI */



  class Vector : public VectorBase
  {
  public:
    typedef ::types::global_dof_index size_type;

    static const bool supports_distributed_data = false;

    Vector () DEAL_II_DEPRECATED;

    explicit Vector (const size_type n) DEAL_II_DEPRECATED;

    explicit Vector (const Epetra_Map &partitioning) DEAL_II_DEPRECATED;

    explicit Vector (const IndexSet &partitioning,
                     const MPI_Comm &communicator = MPI_COMM_WORLD) DEAL_II_DEPRECATED;

    explicit Vector (const VectorBase &V) DEAL_II_DEPRECATED;

    template <typename Number>
    explicit Vector (const ::Vector<Number> &v) DEAL_II_DEPRECATED;

    void reinit (const size_type n,
                 const bool      omit_zeroing_entries = false);

    void reinit (const Epetra_Map &input_map,
                 const bool        omit_zeroing_entries = false);

    void reinit (const IndexSet   &input_map,
                 const MPI_Comm   &communicator = MPI_COMM_WORLD,
                 const bool        omit_zeroing_entries = false);

    void reinit (const VectorBase &V,
                 const bool        omit_zeroing_entries = false,
                 const bool        allow_different_maps = false);

    Vector &operator= (const TrilinosScalar s);

    Vector &operator= (const MPI::Vector &v);

    template <typename Number>
    Vector &operator= (const ::Vector<Number> &v);

    Vector &operator= (const Vector &v);

    void update_ghost_values () const;
  };



// ------------------- inline and template functions --------------


  inline
  void swap (Vector &u, Vector &v)
  {
    u.swap (v);
  }


#ifndef DOXYGEN

  template <typename number>
  Vector::Vector (const ::Vector<number> &v)
  {
    Epetra_LocalMap map ((TrilinosWrappers::types::int_type)v.size(), 0, Utilities::Trilinos::comm_self());
    vector.reset (new Epetra_FEVector(map));
    *this = v;
  }



  inline
  Vector &
  Vector::operator= (const TrilinosScalar s)
  {
    VectorBase::operator= (s);

    return *this;
  }



  template <typename Number>
  Vector &
  Vector::operator= (const ::Vector<Number> &v)
  {
    if (size() != v.size())
      {
        vector.reset();

        Epetra_LocalMap map ((TrilinosWrappers::types::int_type)v.size(), 0,
                             Utilities::Trilinos::comm_self());
        vector.reset (new Epetra_FEVector(map));
      }

    const Epetra_Map &map = vector_partitioner();
    const TrilinosWrappers::types::int_type size = map.NumMyElements();

    Assert (map.MaxLID() == size-1,
            ExcDimensionMismatch(map.MaxLID(), size-1));

    // Need to copy out values, since the
    // deal.II might not use doubles, so
    // that a direct access is not possible.
    for (TrilinosWrappers::types::int_type i=0; i<size; ++i)
      (*vector)[0][i] = v(i);

    return *this;
  }



  inline
  void
  Vector::update_ghost_values () const
  {}


#endif

} /* namespace TrilinosWrappers */

namespace internal
{
  namespace LinearOperator
  {
    template <typename> class ReinitHelper;

    template<>
    class ReinitHelper<TrilinosWrappers::MPI::Vector>
    {
    public:
      template <typename Matrix>
      static
      void reinit_range_vector (const Matrix &matrix,
                                TrilinosWrappers::MPI::Vector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit(matrix.locally_owned_range_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
      }

      template <typename Matrix>
      static
      void reinit_domain_vector(const Matrix &matrix,
                                TrilinosWrappers::MPI::Vector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit(matrix.locally_owned_domain_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
      }
    };

    template<>
    class ReinitHelper<TrilinosWrappers::Vector>
    {
    public:
      template <typename Matrix>
      static
      void reinit_range_vector (const Matrix &matrix,
                                TrilinosWrappers::Vector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit(matrix.locally_owned_range_indices(),
                 matrix.get_mpi_communicator(),
                 omit_zeroing_entries);
      }

      template <typename Matrix>
      static
      void reinit_domain_vector(const Matrix &matrix,
                                TrilinosWrappers::Vector &v,
                                bool omit_zeroing_entries)
      {
        v.reinit(matrix.locally_owned_domain_indices(),
                 matrix.get_mpi_communicator(),
                 omit_zeroing_entries);
      }
    };

  } /* namespace LinearOperator */
} /* namespace internal */


DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_TRILINOS

/*----------------------------   trilinos_vector.h     ---------------------------*/

#endif
/*----------------------------   trilinos_vector.h     ---------------------------*/