ComputeSparsityPattern.cpp File Reference

#include "ComputeSparsityPattern.h"
#include <numeric>
#include <range/v3/algorithm/sort.hpp>
#include <range/v3/algorithm/unique.hpp>
#include <range/v3/range/conversion.hpp>
#include <range/v3/view/transform.hpp>
#include "LocalToGlobalIndexMap.h"
#include "MeshLib/NodeAdjacencyTable.h"
#include "MeshLib/NodePartitionedMesh.h"

Include dependency graph for ComputeSparsityPattern.cpp:

Go to the source code of this file.

Namespaces
namespace	NumLib

Functions
GlobalSparsityPattern	computeSparsityPatternPETSc (NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
GlobalSparsityPattern	NumLib::computeSparsityPattern (LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
	Computes a sparsity pattern for the given inputs.

Function Documentation

computeSparsityPatternPETSc()

GlobalSparsityPattern computeSparsityPatternPETSc	(	NumLib::LocalToGlobalIndexMap const &	dof_table,
		MeshLib::Mesh const &	mesh )

Definition at line 18 of file ComputeSparsityPattern.cpp.

{
    assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(&mesh));
 
    MeshLib::NodeAdjacencyTable const node_adjacency_table(mesh);
 
    // getGlobalIndices() returns raw PETSc global indices:
    //   idx >= 0  — owned by this rank, actual PETSc global row/col index
    //   idx < 0   — ghost node on this rank (owned by another rank);
    //               the actual global index is stored negated
    auto const global_idcs =
        MeshLib::views::meshLocations(mesh, MeshLib::MeshItemType::Node) |
        ranges::views::transform([&](auto&& l)
                                 { return dof_table.getGlobalIndices(l); }) |
        ranges::to<std::vector>();
 
    auto const n_local =
        static_cast<GlobalIndexType>(dof_table.dofSizeWithoutGhosts());
 
    // Each rank owns a contiguous range of global indices; global_start is the
    // first index owned here, obtained by an exclusive prefix-sum of the local
    // sizes. It maps a global row index to the corresponding local row index.
    //
    // MPIU_INT is not an MPI standard datatype but a PETSc macro (defined in
    // petscsys.h) that expands to the MPI datatype matching the width of
    // PetscInt, i.e. MPI_INT or MPI_INT64_T depending on the PETSc
    // configuration. It is used rather than the BaseLib::MPI wrappers because
    // BaseLib::MPI::mpiType() has no mapping for a 64-bit PetscInt and would
    // silently pick the wrong datatype.
    GlobalIndexType global_start = 0;
    MPI_Exscan(&n_local, &global_start, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD);
    // Total global DOF count needed to decode the special ghost encoding where
    // global index 0 is stored as -num_global_dof.
    GlobalIndexType num_global_dof = 0;
    MPI_Allreduce(&n_local, &num_global_dof, 1, MPIU_INT, MPI_SUM,
                  PETSC_COMM_WORLD);
 
    // Collect all (row, col) pairs into per-local-row lists of global column
    // indices; they are sorted and deduplicated below.
    //
    // A row holds only a few dozen columns, so a flat vector plus one
    // sort/unique per row is used instead of std::set: appending is a
    // contiguous write without a per-element node allocation, and the sort
    // is on data that already sits in cache. The result is the same -- sorted
    // and unique -- as std::set gives, which is what the CSR pattern needs.
    std::vector<std::vector<GlobalIndexType>> col_lists(n_local);
 
    // Decode a raw PETSc global index (possibly negative ghost encoding) to the
    // actual non-negative global index.
    auto const decode_col = [&](GlobalIndexType const raw) -> GlobalIndexType
    {
        if (raw >= 0)
        {
            return raw;
        }
        // Ghost encoding: actual index k is stored as -k, except k==0 which is
        // stored as -num_global_dof.
        return (-raw == num_global_dof) ? 0 : -raw;
    };
 
    // Collect a (row, col) pair into the per-row column set.
    //   row >= 0 (non-negative)  — owned by this rank
    //   row < 0                  — ghost row, not owned here → skip
    auto const collect =
        [&](GlobalIndexType const row, GlobalIndexType const col)
    {
        if (row < 0)
        {
            return;
        }
        // A non-negative raw index is owned by this rank, and an owned index
        // lies in this rank's contiguous range [global_start, global_start +
        // n_local), so the local row is always in range here.
        GlobalIndexType const local_row = row - global_start;
        assert(local_row >= 0 && local_row < n_local);
        col_lists[local_row].push_back(decode_col(col));
    };
 
    // Standard mesh-adjacency contribution.
    // Note: getAdjacentNodes() already includes node n itself (the adjacency
    // table is built from all element nodes, including self), so there is no
    // need for a separate self-node loop.
    for (std::size_t n = 0; n < mesh.getNumberOfNodes(); ++n)
    {
        auto const& adj_nodes = node_adjacency_table.getAdjacentNodes(n);
        for (auto const row_idx : global_idcs[n])
        {
            for (auto const adj : adj_nodes)
            {
                for (auto const col_idx : global_idcs[adj])
                {
                    collect(row_idx, col_idx);
                }
            }
        }
    }
 
    // Sort and deduplicate each row, so that col_idx ends up sorted within each
    // row as PETScSparsityPattern documents.
    std::size_t number_of_nonzeros = 0;
    for (auto& cols : col_lists)
    {
        ranges::sort(cols);
        cols.erase(ranges::unique(cols), cols.end());
        number_of_nonzeros += cols.size();
    }
 
    // Build the CSR sparsity pattern from the collected column lists.
    GlobalSparsityPattern sparsity_pattern;
    sparsity_pattern.row_ptr.resize(n_local + 1);
    sparsity_pattern.row_ptr[0] = 0;
    sparsity_pattern.col_idx.reserve(number_of_nonzeros);
 
    for (GlobalIndexType local_row = 0; local_row < n_local; ++local_row)
    {
        auto const& cols = col_lists[local_row];
        sparsity_pattern.col_idx.insert(sparsity_pattern.col_idx.end(),
                                        cols.begin(), cols.end());
        sparsity_pattern.row_ptr[local_row + 1] =
            sparsity_pattern.row_ptr[local_row] +
            static_cast<GlobalIndexType>(cols.size());
    }
 
    return sparsity_pattern;
}

References MathLib::PETScSparsityPattern::col_idx, NumLib::LocalToGlobalIndexMap::dofSizeWithoutGhosts(), MeshLib::NodeAdjacencyTable::getAdjacentNodes(), NumLib::LocalToGlobalIndexMap::getGlobalIndices(), MeshLib::Mesh::getNumberOfNodes(), MeshLib::views::meshLocations(), MeshLib::Node, and MathLib::PETScSparsityPattern::row_ptr.

Referenced by NumLib::computeSparsityPattern().