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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#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/petsc_vector_base.h>
#    include <deal.II/lac/vector_operation.h>

#    include <petscmat.h>

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

DEAL_II_NAMESPACE_OPEN

template <typename Matrix>
class BlockMatrixBase;


namespace PETScWrappers
{
  // forward declarations
  class MatrixBase;

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

        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:
      using size_type = types::global_dof_index;

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

  } // namespace MatrixIterators


  class MatrixBase : public Subscriptor
  {
  public:
    using const_iterator = MatrixIterators::const_iterator;

    using size_type = types::global_dof_index;

    using value_type = PetscScalar;

    MatrixBase();

    MatrixBase(const MatrixBase &) = delete;

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

    virtual ~MatrixBase() override;

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

  } // namespace MatrixIterators



  // 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
} // namespace PETScWrappers


DEAL_II_NAMESPACE_CLOSE


#  endif // DEAL_II_WITH_PETSC

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