petsc_parallel_block_sparse_matrix.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

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

#ifdef DEAL_II_WITH_PETSC

DEAL_II_NAMESPACE_OPEN

namespace PETScWrappers
{
  namespace MPI
  {
    BlockSparseMatrix &
    BlockSparseMatrix::operator=(const BlockSparseMatrix &m)
    {
      BaseClass::operator=(m);

      return *this;
    }


    void
    BlockSparseMatrix::reinit(const size_type n_block_rows,
                              const size_type n_block_columns)
    {
      // first delete previous content of
      // the subobjects array
      clear();

      // then resize. set sizes of blocks to
      // zero. user will later have to call
      // collect_sizes for this
      this->sub_objects.reinit(n_block_rows, n_block_columns);
      this->row_block_indices.reinit(n_block_rows, 0);
      this->column_block_indices.reinit(n_block_columns, 0);

      // and reinitialize the blocks
      for (size_type r = 0; r < this->n_block_rows(); ++r)
        for (size_type c = 0; c < this->n_block_cols(); ++c)
          {
            BlockType *p            = new BlockType();
            this->sub_objects[r][c] = p;
          }
    }

    void
    BlockSparseMatrix::reinit(const std::vector<IndexSet> &      rows,
                              const std::vector<IndexSet> &      cols,
                              const BlockDynamicSparsityPattern &bdsp,
                              const MPI_Comm &                   com)
    {
      Assert(rows.size() == bdsp.n_block_rows(), ExcMessage("invalid size"));
      Assert(cols.size() == bdsp.n_block_cols(), ExcMessage("invalid size"));


      clear();
      this->sub_objects.reinit(bdsp.n_block_rows(), bdsp.n_block_cols());

      std::vector<types::global_dof_index> row_sizes;
      for (unsigned int r = 0; r < bdsp.n_block_rows(); ++r)
        row_sizes.push_back(bdsp.block(r, 0).n_rows());
      this->row_block_indices.reinit(row_sizes);

      std::vector<types::global_dof_index> col_sizes;
      for (unsigned int c = 0; c < bdsp.n_block_cols(); ++c)
        col_sizes.push_back(bdsp.block(0, c).n_cols());
      this->column_block_indices.reinit(col_sizes);

      for (unsigned int r = 0; r < this->n_block_rows(); ++r)
        for (unsigned int c = 0; c < this->n_block_cols(); ++c)
          {
            Assert(rows[r].size() == bdsp.block(r, c).n_rows(),
                   ExcMessage("invalid size"));
            Assert(cols[c].size() == bdsp.block(r, c).n_cols(),
                   ExcMessage("invalid size"));

            BlockType *p = new BlockType();
            p->reinit(rows[r], cols[c], bdsp.block(r, c), com);
            this->sub_objects[r][c] = p;
          }

      collect_sizes();
    }

    void
    BlockSparseMatrix::reinit(const std::vector<IndexSet> &      sizes,
                              const BlockDynamicSparsityPattern &bdsp,
                              const MPI_Comm &                   com)
    {
      reinit(sizes, sizes, bdsp, com);
    }



    void
    BlockSparseMatrix::collect_sizes()
    {
      BaseClass::collect_sizes();
    }

    std::vector<IndexSet>
    BlockSparseMatrix::locally_owned_domain_indices() const
    {
      std::vector<IndexSet> index_sets;

      for (unsigned int i = 0; i < this->n_block_cols(); ++i)
        index_sets.push_back(this->block(0, i).locally_owned_domain_indices());

      return index_sets;
    }

    std::vector<IndexSet>
    BlockSparseMatrix::locally_owned_range_indices() const
    {
      std::vector<IndexSet> index_sets;

      for (unsigned int i = 0; i < this->n_block_rows(); ++i)
        index_sets.push_back(this->block(i, 0).locally_owned_range_indices());

      return index_sets;
    }

    const MPI_Comm &
    BlockSparseMatrix::get_mpi_communicator() const
    {
      return block(0, 0).get_mpi_communicator();
    }

  } // namespace MPI
} // namespace PETScWrappers



DEAL_II_NAMESPACE_CLOSE

#endif