欢迎加入 CSR 存图神教

· · 算法·理论

说在前面

链式前向星访问太慢了?std::vector 模拟邻接表建边太慢了?如果您也面临这样的问题,那么欢迎加入 CSR 存图神教!

CSR 存图法具有 std::vector 模拟邻接表的访问速度和链式前向星的加边速度,同时具有很高的可扩展性,并且实现简单。下面本文将介绍这一存图法的原理和实现,以及一些常见的扩展。最后有一份模板。

原理实现

CSR 存图法的原理非常简单:统计每个点的出度 \mathrm{degree}_u,对它做前缀和,然后将点 u 的出边放进区间 [\mathrm{degree}_{u-1},\mathrm{degree}_u) 中。仿照计数排序,可以写出如下代码:

int U[M], V[M], G[M];
int degree[N + 1];

for (int i = 0; i < m; ++i) {
    cin >> U[i] >> V[i]; // 读入有向边 u -> v
    ++degree[U[i]];      // 统计 u 点出度
}
for (int i = 0; i < n; ++i)
    degree[i + 1] += degree[i]; // 原地前缀和
for (int i = 0; i < m; ++i)
    G[--degree[U[i]]] = V[i];   // 加边

for (int u = 0; u < n; ++u) {
    for (int i = degree[u]; i < degree[u + 1]; ++i) {
        int v = G[i]; // 访问出边 u -> v
        /* do something */
    }
}

以上代码要求点编号 0-indexed 且为有向无权图。注意到在建图过程中,这种方法没有涉及类似 std::vector 的动态扩容;在访问过程中,这种方法没有涉及链式前向星的非连续访问,因此效率比它们都要高很多。

效率测试

:::info[测试环境]{open} 测试参考"如何存图"。在 Luogu IDE 进行测试,版本选择为 C++14 (GCC 9) O2。

编译开关为 g++ -x c++ -std=c++14 -fPIC -DONLINE_JUDGE -Wall -fno-asm -lm -march=native -O2。

洛谷评测机为 Intel(R) Xeon(R) Platinum 8369HC @ 3.30GHz(睿频 3.80GHz)。

测试方法为 5 次结果取平均。单位为毫秒。

:::info[测试代码]

#include <bits/stdc++.h>

constexpr int N = 1000000, M = 2000000;

struct Timer {
    std::chrono::high_resolution_clock::time_point beginTime;
    Timer() : beginTime(std::chrono::high_resolution_clock::now()) {}
    long long operator()() const { return std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now() - beginTime).count(); }
};

struct VectorList {
    std::vector<std::pair<int, int>> graph[N];
    inline void addEdge(int u, int v, int w) { graph[u].emplace_back(v, w); }

    long long build() {
        std::mt19937 rng;
        std::uniform_int_distribution<int> distribute(0, N - 1);
        Timer timer;
        std::for_each(graph, graph + N, [](std::vector<std::pair<int, int>>& i) { i.reserve(32); });
        for (int i = M; i--; addEdge(distribute(rng), distribute(rng), rng()));
        return timer();
    }

    long long traverse() {
        std::mt19937 rng;
        std::vector<int> order(N);
        std::iota(order.begin(), order.end(), 0);
        std::shuffle(order.begin(), order.end(), rng);
        int something = 0;
        Timer timer;
        for (int i : order)
            for (auto edge : graph[i])
                something ^= edge.first ^ edge.second & rng();
        return timer();
    }
};

struct ForwardStar {
    struct Edge { int destination, weight; Edge *next; } edges[M], *head[N], *current = edges - 1;
    inline void addEdge(int u, int v, int w) { *++current = Edge{v, w, head[u]}, head[u] = current; }

    long long build() {
        Timer timer;
        std::mt19937 rng;
        std::uniform_int_distribution<int> distribute(0, N - 1);
        for (int i = M; i--; addEdge(distribute(rng), distribute(rng), rng()));
        return timer();
    }

    long long traverse() {
        std::mt19937 rng;
        std::vector<int> order(N);
        std::iota(order.begin(), order.end(), 0);
        std::shuffle(order.begin(), order.end(), rng);
        int something = 0;
        Timer timer;
        for (int i : order)
            for (auto edge = head[i]; edge; edge = edge->next)
                something ^= edge->destination ^ edge->weight & rng();
        return timer();
    }
};

struct CsrGraph {
    int U[M], V[M], W[M];
    std::pair<int, int> G[M];
    int degree[N + 1];

    long long build() {
        Timer timer;
        std::mt19937 rng;
        std::uniform_int_distribution<int> distribute(0, N - 1);
        for (int i = 0; i < M; ++i) {
            U[i] = distribute(rng), V[i] = distribute(rng), W[i] = rng();
            ++degree[U[i]];
        }
        for (int i = 0; i < N; ++i)
            degree[i + 1] += degree[i];
        for (int i = 0; i < M; ++i)
            G[--degree[U[i]]] = std::pair{V[i], W[i]};
        return timer();
    }

    long long traverse() {
        std::mt19937 rng;
        std::vector<int> order(N);
        std::iota(order.begin(), order.end(), 0);
        std::shuffle(order.begin(), order.end(), rng);
        int something = 0;
        Timer timer;
        for (int i : order)
            for (int o = degree[i]; o < degree[i + 1]; ++o)
                something ^= G[o].first ^ G[o].second & rng();
        return timer();
    }
};

int main() {
    static VectorList vectorList;
    static ForwardStar forwardStar;
    static CsrGraph csrGraph;
    std::cout << "N = " << N << ", M = " << M << '\n';
    std::cout << "vectorList: build time used: " << vectorList.build() << "ms.\n";
    std::cout << "vectorList: traverse time used: " << vectorList.traverse() << "ms.\n";
    std::cout << "forwardStar: build time used: " << forwardStar.build() << "ms.\n";
    std::cout << "forwardStar: traverse time used: " << forwardStar.traverse() << "ms.\n";
    std::cout << "csrGraph: build time used: " << csrGraph.build() << "ms.\n";
    std::cout << "csrGraph: traverse time used: " << csrGraph.traverse() << "ms.\n";
    return 0;
}

:::

稀疏图建边 稀疏图遍历 稠密图建边 稠密图遍历
std::vector 邻接表 326.2 33.8 1386.4 93.8
链式前向星 143.6 98.4 713.4 849.6
CSR 存图法 142.8 18.4 718.8 65.0

几乎完全爆杀另外两个!

扩展使用

只能存有向无权图显然不够,下面是一些扩展使用。

有权情况

额外记录边权数组 W[M] 同时将 G[M] 的类型改为 std::pair 即可。

int U[M], V[M], W[M]; // W[i] 存第 i 条边的边权
pair<int, int> G[M];  // .first: 终点;.second:边权
int degree[N + 1];

for (int i = 0; i < M; ++i) {
    cin >> U[i] >> V[i] >> W[i];
    ++degree[U[i]];
}
for (int i = 0; i < N; ++i)
    degree[i + 1] += degree[i];
for (int i = 0; i < M; ++i)
    G[--degree[U[i]]] = std::make_pair(V[i], W[i]); // 加入边权


for (int u = 0; u < n; ++u) {
    for (int i = degree[u]; i < degree[u + 1]; ++i) {
        int v = G[i].first, w = G[i].second; // 访问出边 u -w-> v
        /* do something */
    }
}

找反向边

常见于网络流算法。双向加边时,记录反向边的指针,存到 G[M] 中。

// 模拟 HLPP 加边
int U[M], V[M], W[M];

struct Edge {
    int v, w;
    Edge* rev; // 指向反向边
} G[M];

int degree[N + 1];

for (int i = 0; i < M; ++i) {
    cin >> U[i] >> V[i] >> W[i];
    ++degree[U[i]], ++degree[V[i]]; // 双向加边,两边都统计出度
}
for (int i = 0; i < N; ++i)
    degree[i + 1] += degree[i];
for (int i = 0; i < M; ++i) {
    Edge* u = G + --degree[U[i]]; // 边 u -> v 的指针
    Edge* v = G + --degree[V[i]]; // 边 v -> u 的指针
    *u = Edge{V[i], W[i], v};
    *v = Edge{U[i], 0, u}; // 网络流反向建边
}


for (int u = 0; u < n; ++u) {
    for (int i = degree[u]; i < degree[u + 1]; ++i) {
        int v = G[i].v, w = G[i].w;             // 访问边 u -> v
        int rv = G[i].rev->v, rw = G[i].rev->w; // 访问边 v -> u
        /* do something */
    }
}

模板封装

下面是一份简单的 CSR 存图模板类,支持像 std::vector 模拟邻接表一样的出边访问方式(即 for (auto [v, w] : graph[u]))。需要 C++20 标准,图为 0-indexed 编号,默认有向,void 为无权图特化。

:::info[code]

#include <ranges>
#include <span>
#include <numeric>
#include <vector>

template <typename Weight>
class CsrGraph {
    std::vector<std::size_t> degree;
    std::vector<std::pair<std::size_t, Weight>> graph;

public:
    template <std::ranges::forward_range EdgeList> requires std::same_as<
        std::remove_cvref_t<std::ranges::range_value_t<EdgeList>>,
        std::tuple<std::size_t, std::size_t, Weight>
    >
    CsrGraph(std::size_t n, std::size_t m, const EdgeList& edgeList) : degree(++n), graph(m) {
        for (auto [u, v, w] : edgeList) ++degree[u];
        std::partial_sum(degree.begin(), degree.end(), degree.data());
        for (auto [u, v, w] : edgeList) graph[--degree[u]] = std::pair{v, w};
    }

    auto operator[](std::size_t i) -> std::span<std::pair<std::size_t, Weight>> {
        return std::span{graph.data() + degree[i], graph.data() + degree[i + 1]};
    }

    auto operator[](std::size_t i) const -> std::span<const std::pair<std::size_t, Weight>> {
        return std::span{graph.data() + degree[i], graph.data() + degree[i + 1]};
    }
};

template <>
class CsrGraph<void> {
    std::vector<std::size_t> degree;
    std::vector<std::size_t> graph;

public:
    template <std::ranges::forward_range EdgeList> requires std::same_as<
        std::remove_cvref_t<std::ranges::range_value_t<EdgeList>>,
        std::tuple<std::size_t, std::size_t>
    >
    CsrGraph(std::size_t n, std::size_t m, const EdgeList& edgeList) : degree(++n), graph(m) {
        for (auto [u, v] : edgeList) ++degree[u];
        std::partial_sum(degree.begin(), degree.end(), degree.data());
        for (auto [u, v] : edgeList) graph[--degree[u]] = v;
    }

    auto operator[](std::size_t i) -> std::span<std::size_t> {
        return std::span{graph.data() + degree[i], graph.data() + degree[i + 1]};
    }

    auto operator[](std::size_t i) const -> std::span<const std::size_t> {
        return std::span{graph.data() + degree[i], graph.data() + degree[i + 1]};
    }
};

:::

参考资料