petsc_matrix_base.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2018 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_petsc_matrix_base_h
#define dealii_petsc_matrix_base_h


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

#ifdef DEAL_II_WITH_PETSC

#  include <deal.II/base/subscriptor.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/petsc_compatibility.h>
#  include <deal.II/lac/vector_operation.h>
#  include <deal.II/lac/petsc_vector_base.h>

#  include <petscmat.h>

#  include <vector>
#  include <cmath>
#  include <memory>

DEAL_II_NAMESPACE_OPEN

template <typename Matrix> class BlockMatrixBase;


namespace PETScWrappers
{
  // forward declarations
  class MatrixBase;

  namespace MatrixIterators
  {
    class const_iterator
    {
    private:
      class Accessor
      {
      public:
        typedef types::global_dof_index size_type;

        Accessor (const MatrixBase *matrix,
                  const size_type   row,
                  const size_type   index);

        size_type row() const;

        size_type index() const;

        size_type column() const;

        PetscScalar value() const;

        DeclException0 (ExcBeyondEndOfMatrix);
        DeclException3 (ExcAccessToNonlocalRow,
                        int, int, int,
                        << "You tried to access row " << arg1
                        << " of a distributed matrix, but only rows "
                        << arg2 << " through " << arg3
                        << " are stored locally and can be accessed.");

      private:
        mutable MatrixBase *matrix;

        size_type a_row;

        size_type a_index;

        std::shared_ptr<const std::vector<size_type> > colnum_cache;

        std::shared_ptr<const std::vector<PetscScalar> > value_cache;

        void visit_present_row ();

        friend class const_iterator;
      };

    public:
      typedef types::global_dof_index size_type;

      const_iterator (const MatrixBase *matrix,
                      const size_type   row,
                      const size_type   index);

      const_iterator &operator++ ();

      const_iterator operator++ (int);

      const Accessor &operator* () const;

      const Accessor *operator-> () const;

      bool operator == (const const_iterator &) const;
      bool operator != (const const_iterator &) const;

      bool operator < (const const_iterator &) const;

      DeclException2 (ExcInvalidIndexWithinRow,
                      int, int,
                      << "Attempt to access element " << arg2
                      << " of row " << arg1
                      << " which doesn't have that many elements.");

    private:
      Accessor accessor;
    };

  }


  class MatrixBase : public Subscriptor
  {
  public:
    typedef MatrixIterators::const_iterator const_iterator;

    typedef types::global_dof_index size_type;

    typedef PetscScalar value_type;

    MatrixBase ();

    MatrixBase(const MatrixBase &) = delete;

    MatrixBase &operator=(const MatrixBase &) = delete;

    virtual ~MatrixBase ();

    MatrixBase &
    operator = (const value_type d);
    void clear ();

    void set (const size_type   i,
              const size_type   j,
              const PetscScalar value);

    void set (const std::vector<size_type>  &indices,
              const FullMatrix<PetscScalar> &full_matrix,
              const bool                     elide_zero_values = false);

    void set (const std::vector<size_type>  &row_indices,
              const std::vector<size_type>  &col_indices,
              const FullMatrix<PetscScalar> &full_matrix,
              const bool                     elide_zero_values = false);

    void set (const size_type                 row,
              const std::vector<size_type >  &col_indices,
              const std::vector<PetscScalar>  &values,
              const bool                      elide_zero_values = false);

    void set (const size_type    row,
              const size_type    n_cols,
              const size_type   *col_indices,
              const PetscScalar  *values,
              const bool         elide_zero_values = false);

    void add (const size_type   i,
              const size_type   j,
              const PetscScalar value);

    void add (const std::vector<size_type>  &indices,
              const FullMatrix<PetscScalar> &full_matrix,
              const bool                     elide_zero_values = true);

    void add (const std::vector<size_type>  &row_indices,
              const std::vector<size_type>  &col_indices,
              const FullMatrix<PetscScalar> &full_matrix,
              const bool                     elide_zero_values = true);

    void add (const size_type                 row,
              const std::vector<size_type>   &col_indices,
              const std::vector<PetscScalar>  &values,
              const bool                      elide_zero_values = true);

    void add (const size_type    row,
              const size_type    n_cols,
              const size_type   *col_indices,
              const PetscScalar  *values,
              const bool         elide_zero_values = true,
              const bool         col_indices_are_sorted = false);

    void clear_row (const size_type   row,
                    const PetscScalar new_diag_value = 0);

    void clear_rows (const std::vector<size_type> &rows,
                     const PetscScalar             new_diag_value = 0);

    void compress (const VectorOperation::values operation);

    PetscScalar operator () (const size_type i,
                             const size_type j) const;

    PetscScalar el (const size_type i,
                    const size_type j) const;

    PetscScalar diag_element (const size_type i) const;

    size_type m () const;

    size_type n () const;

    size_type local_size () const;

    std::pair<size_type, size_type>
    local_range () const;

    bool in_local_range (const size_type index) const;

    virtual const MPI_Comm &get_mpi_communicator () const = 0;

    size_type n_nonzero_elements () const;

    size_type row_length (const size_type row) const;

    PetscReal l1_norm () const;

    PetscReal linfty_norm () const;

    PetscReal frobenius_norm () const;


    PetscScalar matrix_norm_square (const VectorBase &v) const;


    PetscScalar matrix_scalar_product (const VectorBase &u,
                                       const VectorBase &v) const;

    PetscScalar trace () const;

    MatrixBase &operator *= (const PetscScalar factor);

    MatrixBase &operator /= (const PetscScalar factor);


    MatrixBase &add (const PetscScalar factor,
                     const MatrixBase &other);


    DEAL_II_DEPRECATED
    MatrixBase &add (const MatrixBase &other,
                     const PetscScalar factor);

    void vmult (VectorBase       &dst,
                const VectorBase &src) const;

    void Tvmult (VectorBase       &dst,
                 const VectorBase &src) const;

    void vmult_add (VectorBase       &dst,
                    const VectorBase &src) const;

    void Tvmult_add (VectorBase       &dst,
                     const VectorBase &src) const;

    PetscScalar residual (VectorBase       &dst,
                          const VectorBase &x,
                          const VectorBase &b) const;

    const_iterator begin () const;

    const_iterator end () const;

    const_iterator begin (const size_type r) const;

    const_iterator end (const size_type r) const;

    operator Mat () const;

    Mat &petsc_matrix ();

    void transpose ();

    PetscBool
    is_symmetric (const double tolerance = 1.e-12);

    PetscBool
    is_hermitian (const double tolerance = 1.e-12);

    void write_ascii (const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);

    void print (std::ostream &out,
                const bool    alternative_output = false) const;

    std::size_t memory_consumption() const;

    DeclExceptionMsg (ExcSourceEqualsDestination,
                      "You are attempting an operation on two matrices that "
                      "are the same object, but the operation requires that the "
                      "two objects are in fact different.");

    DeclException2 (ExcWrongMode,
                    int, int,
                    << "You tried to do a "
                    << (arg1 == 1 ?
                        "'set'" :
                        (arg1 == 2 ?
                         "'add'" : "???"))
                    << " operation but the matrix is currently in "
                    << (arg2 == 1 ?
                        "'set'" :
                        (arg2 == 2 ?
                         "'add'" : "???"))
                    << " mode. You first have to call 'compress()'.");

  protected:
    Mat matrix;

    VectorOperation::values last_action;

    void prepare_action(const VectorOperation::values new_action);

    void assert_is_compressed();

    void prepare_add();
    void prepare_set();

    void mmult (MatrixBase       &C,
                const MatrixBase &B,
                const VectorBase &V) const;

    void Tmmult (MatrixBase       &C,
                 const MatrixBase &B,
                 const VectorBase &V) const;

  private:

    mutable std::vector<PetscInt> column_indices;

    mutable std::vector<PetscScalar> column_values;


    template <class> friend class ::BlockMatrixBase;
  };



#ifndef DOXYGEN
// -------------------------- inline and template functions ----------------------


  namespace MatrixIterators
  {

    inline
    const_iterator::Accessor::
    Accessor (const MatrixBase *matrix,
              const size_type   row,
              const size_type   index)
      :
      matrix(const_cast<MatrixBase *>(matrix)),
      a_row(row),
      a_index(index)
    {
      visit_present_row ();
    }



    inline
    const_iterator::Accessor::size_type
    const_iterator::Accessor::row() const
    {
      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return a_row;
    }


    inline
    const_iterator::Accessor::size_type
    const_iterator::Accessor::column() const
    {
      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return (*colnum_cache)[a_index];
    }


    inline
    const_iterator::Accessor::size_type
    const_iterator::Accessor::index() const
    {
      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return a_index;
    }


    inline
    PetscScalar
    const_iterator::Accessor::value() const
    {
      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return (*value_cache)[a_index];
    }


    inline
    const_iterator::
    const_iterator(const MatrixBase *matrix,
                   const size_type   row,
                   const size_type   index)
      :
      accessor(matrix, row, index)
    {}



    inline
    const_iterator &
    const_iterator::operator++ ()
    {
      Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());

      ++accessor.a_index;

      // if at end of line: do one step, then cycle until we find a
      // row with a nonzero number of entries
      if (accessor.a_index >= accessor.colnum_cache->size())
        {
          accessor.a_index = 0;
          ++accessor.a_row;

          while ((accessor.a_row < accessor.matrix->m())
                 &&
                 (accessor.a_row < accessor.matrix->local_range().second)
                 &&
                 (accessor.matrix->row_length(accessor.a_row) == 0))
            ++accessor.a_row;

          accessor.visit_present_row();
        }
      return *this;
    }


    inline
    const_iterator
    const_iterator::operator++ (int)
    {
      const const_iterator old_state = *this;
      ++(*this);
      return old_state;
    }


    inline
    const const_iterator::Accessor &
    const_iterator::operator* () const
    {
      return accessor;
    }


    inline
    const const_iterator::Accessor *
    const_iterator::operator-> () const
    {
      return &accessor;
    }


    inline
    bool
    const_iterator::
    operator == (const const_iterator &other) const
    {
      return (accessor.a_row == other.accessor.a_row &&
              accessor.a_index == other.accessor.a_index);
    }


    inline
    bool
    const_iterator::
    operator != (const const_iterator &other) const
    {
      return ! (*this == other);
    }


    inline
    bool
    const_iterator::
    operator < (const const_iterator &other) const
    {
      return (accessor.row() < other.accessor.row() ||
              (accessor.row() == other.accessor.row() &&
               accessor.index() < other.accessor.index()));
    }

  }



  // Inline the set() and add()
  // functions, since they will be
  // called frequently, and the
  // compiler can optimize away
  // some unnecessary loops when
  // the sizes are given at
  // compile time.
  inline
  void
  MatrixBase::set (const size_type   i,
                   const size_type   j,
                   const PetscScalar value)
  {
    AssertIsFinite(value);

    set (i, 1, &j, &value, false);
  }



  inline
  void
  MatrixBase::set (const std::vector<size_type>  &indices,
                   const FullMatrix<PetscScalar> &values,
                   const bool                     elide_zero_values)
  {
    Assert (indices.size() == values.m(),
            ExcDimensionMismatch(indices.size(), values.m()));
    Assert (values.m() == values.n(), ExcNotQuadratic());

    for (size_type i=0; i<indices.size(); ++i)
      set (indices[i], indices.size(), indices.data(), &values(i,0),
           elide_zero_values);
  }



  inline
  void
  MatrixBase::set (const std::vector<size_type>  &row_indices,
                   const std::vector<size_type>  &col_indices,
                   const FullMatrix<PetscScalar> &values,
                   const bool                     elide_zero_values)
  {
    Assert (row_indices.size() == values.m(),
            ExcDimensionMismatch(row_indices.size(), values.m()));
    Assert (col_indices.size() == values.n(),
            ExcDimensionMismatch(col_indices.size(), values.n()));

    for (size_type i=0; i<row_indices.size(); ++i)
      set (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
           elide_zero_values);
  }



  inline
  void
  MatrixBase::set (const size_type                 row,
                   const std::vector<size_type>   &col_indices,
                   const std::vector<PetscScalar>  &values,
                   const bool                      elide_zero_values)
  {
    Assert (col_indices.size() == values.size(),
            ExcDimensionMismatch(col_indices.size(), values.size()));

    set (row, col_indices.size(), col_indices.data(), values.data(),
         elide_zero_values);
  }



  inline
  void
  MatrixBase::set (const size_type    row,
                   const size_type    n_cols,
                   const size_type   *col_indices,
                   const PetscScalar  *values,
                   const bool         elide_zero_values)
  {
    prepare_action(VectorOperation::insert);

    const PetscInt petsc_i = row;
    PetscInt const *col_index_ptr;

    PetscScalar const *col_value_ptr;
    int n_columns;

    // If we don't elide zeros, the pointers are already available...
    if (elide_zero_values == false)
      {
        col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
        col_value_ptr = values;
        n_columns = n_cols;
      }
    else
      {
        // Otherwise, extract nonzero values in each row and get the
        // respective index.
        if (column_indices.size() < n_cols)
          {
            column_indices.resize(n_cols);
            column_values.resize(n_cols);
          }

        n_columns = 0;
        for (size_type j=0; j<n_cols; ++j)
          {
            const PetscScalar value = values[j];
            AssertIsFinite(value);
            if (value != PetscScalar())
              {
                column_indices[n_columns] = col_indices[j];
                column_values[n_columns] = value;
                n_columns++;
              }
          }
        Assert(n_columns <= (int)n_cols, ExcInternalError());

        col_index_ptr = column_indices.data();
        col_value_ptr = column_values.data();
      }

    const PetscErrorCode ierr = MatSetValues (matrix, 1, &petsc_i, n_columns,
                                              col_index_ptr,
                                              col_value_ptr, INSERT_VALUES);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }



  inline
  void
  MatrixBase::add (const size_type   i,
                   const size_type   j,
                   const PetscScalar value)
  {

    AssertIsFinite(value);

    if (value == PetscScalar())
      {
        // we have to check after using Insert/Add in any case to be
        // consistent with the MPI communication model, but we can save
        // some work if the addend is zero. However, these actions are done
        // in case we pass on to the other function.
        prepare_action(VectorOperation::add);

        return;
      }
    else
      add (i, 1, &j, &value, false);
  }



  inline
  void
  MatrixBase::add (const std::vector<size_type>  &indices,
                   const FullMatrix<PetscScalar> &values,
                   const bool                     elide_zero_values)
  {
    Assert (indices.size() == values.m(),
            ExcDimensionMismatch(indices.size(), values.m()));
    Assert (values.m() == values.n(), ExcNotQuadratic());

    for (size_type i=0; i<indices.size(); ++i)
      add (indices[i], indices.size(), indices.data(), &values(i,0),
           elide_zero_values);
  }



  inline
  void
  MatrixBase::add (const std::vector<size_type>  &row_indices,
                   const std::vector<size_type>  &col_indices,
                   const FullMatrix<PetscScalar> &values,
                   const bool                     elide_zero_values)
  {
    Assert (row_indices.size() == values.m(),
            ExcDimensionMismatch(row_indices.size(), values.m()));
    Assert (col_indices.size() == values.n(),
            ExcDimensionMismatch(col_indices.size(), values.n()));

    for (size_type i=0; i<row_indices.size(); ++i)
      add (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
           elide_zero_values);
  }



  inline
  void
  MatrixBase::add (const size_type                 row,
                   const std::vector<size_type>   &col_indices,
                   const std::vector<PetscScalar>  &values,
                   const bool                      elide_zero_values)
  {
    Assert (col_indices.size() == values.size(),
            ExcDimensionMismatch(col_indices.size(), values.size()));

    add (row, col_indices.size(), col_indices.data(), values.data(),
         elide_zero_values);
  }



  inline
  void
  MatrixBase::add (const size_type    row,
                   const size_type    n_cols,
                   const size_type   *col_indices,
                   const PetscScalar *values,
                   const bool         elide_zero_values,
                   const bool          /*col_indices_are_sorted*/)
  {
    (void)elide_zero_values;

    prepare_action(VectorOperation::add);

    const PetscInt petsc_i = row;
    PetscInt const *col_index_ptr;

    PetscScalar const *col_value_ptr;
    int n_columns;

    // If we don't elide zeros, the pointers are already available...
    if (elide_zero_values == false)
      {
        col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
        col_value_ptr = values;
        n_columns = n_cols;
      }
    else
      {
        // Otherwise, extract nonzero values in each row and get the
        // respective index.
        if (column_indices.size() < n_cols)
          {
            column_indices.resize(n_cols);
            column_values.resize(n_cols);
          }

        n_columns = 0;
        for (size_type j=0; j<n_cols; ++j)
          {
            const PetscScalar value = values[j];
            AssertIsFinite(value);
            if (value != PetscScalar())
              {
                column_indices[n_columns] = col_indices[j];
                column_values[n_columns] = value;
                n_columns++;
              }
          }
        Assert(n_columns <= (int)n_cols, ExcInternalError());

        col_index_ptr = column_indices.data();
        col_value_ptr = column_values.data();
      }

    const PetscErrorCode ierr = MatSetValues (matrix, 1, &petsc_i, n_columns,
                                              col_index_ptr, col_value_ptr,
                                              ADD_VALUES);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }






  inline
  PetscScalar
  MatrixBase::operator() (const size_type i,
                          const size_type j) const
  {
    return el(i,j);
  }



  inline
  MatrixBase::const_iterator
  MatrixBase::begin() const
  {
    return const_iterator(this, 0, 0);
  }


  inline
  MatrixBase::const_iterator
  MatrixBase::end() const
  {
    return const_iterator(this, m(), 0);
  }


  inline
  MatrixBase::const_iterator
  MatrixBase::begin(const size_type r) const
  {
    Assert (in_local_range(r),
            ExcIndexRange(r, local_range().first, local_range().second));

    if (row_length(r) > 0)
      return const_iterator(this, r, 0);
    else
      return end (r);
  }


  inline
  MatrixBase::const_iterator
  MatrixBase::end(const size_type r) const
  {
    Assert (in_local_range(r),
            ExcIndexRange(r, local_range().first, local_range().second));

    // place the iterator on the first entry past this line, or at the
    // end of the matrix
    //
    // in the parallel case, we need to put it on the first entry of
    // the first row after the locally owned range. this of course
    // doesn't exist, but we can nevertheless create such an
    // iterator. we need to check whether 'i' is past the locally
    // owned range of rows first, before we ask for the length of the
    // row since the latter query leads to an exception in case we ask
    // for a row that is not locally owned
    for (size_type i=r+1; i<m(); ++i)
      if (i == local_range().second || (row_length(i) > 0))
        return const_iterator(this, i, 0);

    // if there is no such line, then take the
    // end iterator of the matrix
    return end();
  }



  inline
  bool
  MatrixBase::in_local_range (const size_type index) const
  {
    PetscInt begin, end;

    const PetscErrorCode ierr = MatGetOwnershipRange (static_cast<const Mat &>(matrix),
                                                      &begin, &end);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    return ((index >= static_cast<size_type>(begin)) &&
            (index < static_cast<size_type>(end)));
  }



  inline
  void
  MatrixBase::prepare_action(const VectorOperation::values new_action)
  {
    if (last_action == VectorOperation::unknown)
      last_action = new_action;

    Assert (last_action == new_action, ExcWrongMode (last_action, new_action));
  }



  inline
  void
  MatrixBase::assert_is_compressed ()
  {
    // compress() sets the last action to none, which allows us to check if there
    // are pending add/insert operations:
    AssertThrow (last_action == VectorOperation::unknown,
                 ExcMessage("Error: missing compress() call."));
  }



  inline
  void
  MatrixBase::prepare_add()
  {
    prepare_action(VectorOperation::add);
  }



  inline
  void
  MatrixBase::prepare_set()
  {
    prepare_action(VectorOperation::insert);
  }

#endif // DOXYGEN
}


DEAL_II_NAMESPACE_CLOSE


#endif // DEAL_II_WITH_PETSC


/*----------------------------   petsc_matrix_base.h     ---------------------------*/

#endif
/*----------------------------   petsc_matrix_base.h     ---------------------------*/