SparsityPattern.h

// Copyright (C) 2007-2026 Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include <cstddef>
#include <cstdint>
#include <dolfinx/common/MPI.h>
#include <memory>
#include <span>
#include <utility>
#include <vector>
 
namespace dolfinx::common
{
class IndexMap;
}
 
namespace dolfinx::la
{
class SparsityPattern
{
public:
  SparsityPattern(MPI_Comm comm,
                  std::array<std::shared_ptr<const common::IndexMap>, 2> maps,
                  std::array<int, 2> bs);

  SparsityPattern(
      MPI_Comm comm,
      const std::vector<std::vector<const SparsityPattern*>>& patterns,
      const std::array<
          std::vector<
              std::pair<std::reference_wrapper<const common::IndexMap>, int>>,
          2>& maps,
      const std::array<std::vector<int>, 2>& bs);
 
  SparsityPattern(const SparsityPattern& pattern) = delete;

  SparsityPattern(SparsityPattern&& pattern) = default;

  ~SparsityPattern() = default;

  SparsityPattern& operator=(SparsityPattern&& pattern) = default;

  void reserve_blocks(std::size_t num_blocks, std::size_t num_rows,
                      std::size_t num_cols);

  void insert(std::int32_t row, std::int32_t col);

  void insert(std::span<const std::int32_t> rows,
              std::span<const std::int32_t> cols);

  void insert_diagonal(std::span<const std::int32_t> rows);

  void finalize();

  std::shared_ptr<const common::IndexMap> index_map(int dim) const;

  std::vector<std::int64_t> column_indices() const;

  int block_size(int dim) const;

  std::int64_t num_nonzeros() const;

  std::int32_t nnz_diag(std::int32_t row) const;

  std::int32_t nnz_off_diag(std::int32_t row) const;

  std::pair<std::span<const std::int32_t>, std::span<const std::int64_t>>
  graph() const;

  std::span<const std::int32_t> off_diagonal_offsets() const;

  MPI_Comm comm() const;
 
private:
  // MPI communicator
  dolfinx::MPI::Comm _comm;
 
  // Index maps for each dimension
  std::array<std::shared_ptr<const common::IndexMap>, 2> _index_maps;
 
  // Block size
  std::array<int, 2> _bs;
 
  // Non-zero ghost columns in owned rows
  std::vector<std::int64_t> _col_ghosts;
 
  // Owning process of ghost columns in owned rows
  std::vector<std::int32_t> _col_ghost_owners;
 
  // Cache of unassembled entries on owned and unowned (ghost) rows,
  // held until finalize().
  //
  // insert(rows, cols) inserts the outer product of `rows` and `cols`,
  // so storing one (row, column) pair per entry repeats each index
  // rows.size() or cols.size() times over. Blocks are instead kept in
  // the form they arrive in and expanded only in finalize(): for a P1
  // tetrahedron that is 4 + 4 indices per cell rather than 16 + 16.
  std::vector<std::int32_t> _cache_brows, _cache_bcols;
  std::vector<std::int64_t> _cache_boffs_r{0}, _cache_boffs_c{0};
 
  // Cache of square blocks, i.e. those whose row and column index lists
  // are the same span, which is the case whenever the test and trial
  // dofmaps and the cells indexing them coincide. Only the row list is
  // stored, halving both the copy in insert() and the cached bytes.
  //
  // _cache_sbs is the width shared by every cached square block, 0 if
  // none are cached and -1 if they differ. Cell integrals give blocks
  // of a single width, which are indexed by stride, so _cache_soffs
  // stays empty until a block of a different width arrives.
  std::vector<std::int32_t> _cache_srows;
  std::int32_t _cache_sbs = 0;
  std::vector<std::int64_t> _cache_soffs;
 
  // Cache of individually inserted (row, column) pairs (row-major COO)
  std::vector<std::int32_t> _cache_rows;
  std::vector<std::int32_t> _cache_cols;
 
  // Rows with a cached diagonal entry, from insert_diagonal
  std::vector<std::int32_t> _cache_diag;
 
  // Pending reserve_blocks request, as {number of blocks, total row
  // indices, total column indices}. Which of the two block caches the
  // blocks land in is only known once insert() sees the spans, and
  // reserving both maps around 2.5x what is used, so the request is
  // held here until then.
  std::array<std::size_t, 3> _reserve{0, 0, 0};

  void expand_square_offsets();

  std::pair<std::vector<std::int64_t>, std::vector<std::int32_t>>
  bucket_cache(std::int32_t num_rows, std::int32_t num_cols) const;
 
  // Sparsity pattern adjacency data (computed once pattern is
  // finalised). _edges holds the edges (connected dofs). The edges for
  // node i are in the range [_offsets[i], _offsets[i + 1]).
  std::vector<std::int32_t> _edges;
  std::vector<std::int64_t> _offsets;
 
  // Start of off-diagonal (unowned columns) on each row (row-wise)
  std::vector<std::int32_t> _off_diagonal_offsets;
};
} // namespace dolfinx::la