预备知识：Lambda 表达式

语法

Lambda 表达式定义了一个匿名函数，并且可以捕获一定范围内的变量。lambda 表达式的语法形式可简单归纳如下：

1
2
3
4
5
6
// capture 控制外部变量捕获方式。
// params 是参数表。
// opt 是可选修饰，例如 mutable。
// ret 是返回类型。
// body 是函数体。
[capture](params)opt->ret {body;};

其中 capture 是捕获列表，控制表达式可以捕获指定范围的变量；params 是参数表，传入函数要使用的参数；opt 是函数选项，如mutable可以让表达式可以对创建的捕获变量的副本进行修改和成员函数调用；ret 是返回值类型；body 是函数体。

[capture] 捕获范围

通过不同的 capture 配置，我们可以让 Lambda 函数捕获不同范围内的变量：

  • [] 不捕获任何变量。
  • [&] 捕获外部作用域中所有变量，并作为引用在函数体中使用（按引用捕获）。
  • [=] 捕获外部作用域中所有变量，并作为副本在函数体中使用（按值捕获）。
  • [=，&foo] 按值捕获外部作用域中所有变量，并按引用捕获 foo 变量。
  • [bar] 按值捕获 bar 变量，同时不捕获其他变量。
  • [this] 捕获当前类中的 this 指针，让 lambda 表达式拥有和当前类成员函数同样的访问权限。如果已经使用了 & 或者 =，就默认添加此选项。捕获 this 的目的是可以在 lambda 中使用当前类的成员函数和成员变量。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
// 示例类：演示 lambda 捕获类成员和局部变量的规则。
class A {
public:
    // 类成员变量，可通过 this 指针访问。
    int i_ = 0;
    // x、y 是成员函数的局部参数。
    void func(int x, int y) {
        // 错误：空捕获列表不能访问成员变量 i_。
        auto x1 = { return i_; }; // 编译错误，未捕获外部变量
        // 合法：[=] 默认按值捕获局部变量，并隐式捕获 this。
        auto x2 = [=] { return i_ + x + y; }; // 合法，按值捕获所有变量
        // 合法：[&] 默认按引用捕获局部变量，并隐式捕获 this。
        auto x3 = [&] { return i_ + x + y; }; // 合法，按引用捕获所有变量
        // 合法：显式捕获 this 后可访问成员变量 i_。
        auto x4 = [this] { return i_; }; // 合法，捕获当前对象指针
        // 错误：只捕获 this，不能访问未捕获的局部变量 x 和 y。
        auto x5 = [this] { return i_ + x + y; }; // 编译错误，未捕获局部变量x和y
        // 合法：捕获 this 访问成员变量，按值捕获 x 和 y。
        auto x6 = [this, x, y] { return i_ + x + y; }; // 合法，精确指定捕获对象
        // 合法：通过 this 修改当前对象的成员变量。
        auto x7 = [this] { return i_++; }; // 合法，通过对象指针修改成员变量
    }
};

// 外部局部变量。
int a = 0, b = 1;
// 错误：空捕获列表不能访问 a。
auto f1 = { return a; }; // 编译错误，未捕获外部变量
// 合法：按引用捕获 a，可修改原变量。
auto f2 = [&] { return a++; }; // 合法，按引用捕获并修改
// 合法：按值捕获 a，只读取副本。
auto f3 = [=] { return a; }; // 合法，按值捕获读取
// 错误：按值捕获的副本默认只读，不能执行自增。
auto f4 = [=] { return a++; }; // 编译错误，按值捕获的副本默认不可变
// 错误：只捕获 a，不能访问 b。
auto f5 = [a] { return a + b; }; // 编译错误，局部未捕获变量b
// 合法：a 按值捕获，b 按引用捕获并可修改。
auto f6 = [a, &b] { return a + (b++); }; // 合法，精确控制读写权限
// 合法：默认按值捕获其余变量，但 b 按引用捕获。
auto f7 = [=, &b] { return a + (b++); }; // 合法，混合捕获模式

应用

以统计集合中的偶数数量为例，传统的方法需要定义包含状态变量的结构体，重载函数调用运算符，最后在算法接口中进行实例化调用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
class CountEven {
    // 保存外部计数器引用，使函数对象能够修改调用方变量。
    int& count;
public:
    // 构造时绑定外部 even_count。
    CountEven(int& count) : count(count) {}
    // for_each 每访问一个元素，就调用一次 operator()。
    void operator() (int val) {
        // val & 1 为 0 表示偶数。
        if (!(val & 1)) {
            // 通过引用修改外部计数器。
            ++count;
        }
    }
};

// 待统计的输入序列。
std::vector<int> v = {1, 2, 3, 4, 5, 6};
// 统计结果变量，由 CountEven 持有引用。
int even_count = 0;
// 将函数对象传给标准算法，遍历区间 [begin, end)。
for_each(v.begin(), v.end(), CountEven(even_count));
std::cout << "The number of even is " << even_count << std::endl;

采用 Lambda 后，逻辑被集中在调用处，通过按引用捕获外部计数器，减少了样板代码。这种方式在并行计算逻辑分发中，能够提高程序的并发组织效率与可维护性：

1
2
3
4
5
6
7
8
9
10
11
12
13
// 待统计的输入序列。
std::vector<int> v = {1, 2, 3, 4, 5, 6};
// 统计结果变量。
int even_count = 0;
// 按引用捕获 even_count，使 lambda 内部可以累加外部计数器。
for_each(v.begin(), v.end(), [&even_count] (int val) {
    // 对每个元素执行一次判断，偶数则累加。
    if (!(val & 1)) {
        ++even_count;
    }
});
// 输出最终偶数数量。
std::cout << "The number of even is " << even_count << std::endl;

PageRank 算法

PageRank

互联网上存在大量互相依赖或重复的网页，且规模会随着时间持续增长。依靠纯文本关键词匹配的搜索算法，容易受到低质量信息以及恶意伪造权重网页的干扰，搜索结果的精确度与信噪比难以满足需求。通过分析超链接构建的有向图拓扑特征，利用节点之间的连接关系来计算客观的权重指标，成为提升信息检索质量的突破口。

因此要找到一个量来定义一个页面的“影响力”或者“重要性”。要评估一个页面的重要性，指向这个页面的页面数量就是一个很直观的评估指标——依赖这个页面的页面数量越多，这个页面的重要性就越大。然而还有一个问题，由于不同页面的重要性是不一致的，我们不能用被指向的数量来简单决定页面的重要性，这样一个被 10 个重要页面指向的页面，其重要性甚至不如被 100 个垃圾页面指向的页面。解决方法很简单，计算的时候带上跳转前页面的权重，再乘一个进入被指向页面的概率即可，因此有了计算 PageRank 的公式：

PR(u)=∑v∈D(u)PR(v)∣S(v)∣PR(u)=\sum_{v \in D(u)} \frac{PR(v)}{|S(v)|}

其中 PR(u)PR(u) 代表网页 uu 的 rank 值；PR(v)PR(v) 代表网页 vv 的 rank 值；D(u)D(u) 代表指向 uu 的所有网页，也就是网页图中 uu 的入边源点集合；S(v)S(v) 代表从 vv 出发指向的网页集合，∣S(v)∣|S(v)| 就是该集合的大小，其倒数代表进入被指向页面的概率。

另外一个要考虑的问题是用户在访问某个页面时，会不会通过这个页面进入另一个页面并访问。这个指标直接决定了计算 PageRank 时考虑的网页之间的连接，是否真的“有效”——如果完全没有用户从这个路径去访问另一个页面，那么计算 PageRank 时考虑这些连接就是多余的。因此再引入一个阻尼系数 dd，代表用户从一个页面点击通往另一个网页的链接的概率；相应的，1−d1-d 就是停止点击的概率，那么公式变成：

PR(u)=(1−d)+d×∑v∈D(u)PR(v)∣S(v)∣PR(u)=(1-d) + d \times \sum_{v \in D(u)} \frac{PR(v)}{|S(v)|}

那么，用户点击链接在网页间跳转的概率较大时，就赋予 ∑v∈D(u)PR(v)∣S(v)∣\sum_{v \in D(u)} \frac{PR(v)}{|S(v)|} 较大的权重，让网页间跳转连接的分布来控制网页的重要性；如果用户不倾向于在网页间跳转，就用相同的 1−d1-d 作为网页重要性的“兜底”——权重趋于降低时，不同网页的重要性趋于一致。

用一个简单的示例来演示 PageRank 的计算：

初始状态下，系统将总权重平均分配给所有节点，因此v1到v4的初始权重值均为 0.25，系统根据图的入边关系和源点的出度进行第一轮计算：

  • v1 的计算：v1 接收来自某出度为1的节点（传递全部权重）以及某出度为2的节点（传递一半权重）的值。计算过程为：PR(v1)=0.25×1+0.25×12=0.375PR(v1)=0.25\times 1+0.25\times \frac{1}{2}=0.375（PPT中四舍五入为0.37）。
  • v2 的计算：v2 仅接收来自某出度为3的节点的三分之一权重。计算过程为：PR(v2)=0.25×13≈0.08PR(v2)=0.25\times \frac{1}{3}\approx 0.08。
  • v3 的计算：v3 接收来自三个不同节点的权重传递（源节点出度分别为3、2、2）。计算过程为：PR(v3)=0.25×13+0.25×12+0.25×12≈0.33PR(v3)=0.25\times \frac{1}{3}+0.25\times \frac{1}{2}+0.25\times \frac{1}{2}\approx 0.33。
  • v4 的计算：v4 接收来自两个不同节点的权重传递（源节点出度分别为3、2）。计算过程为：PR(v4)=0.25×13+0.25×12≈0.20PR(v4)=0.25\times \frac{1}{3}+0.25\times \frac{1}{2}\approx 0.20。

完成首轮更新后，将产生的新权重组（0.37, 0.08, 0.33, 0.20）代回公式，触发第二次迭代，并依此类推。随着迭代步数增加，数据呈现如下演化趋势 ：

  • 第2次迭代：v1为0.43，v2为0.12，v3为0.27，v4为0.16。
  • 第3次迭代：v1为0.35，v2为0.14，v3为0.29，v4为0.20。
  • 第4次迭代：v1为0.39，v2为0.11，v3为0.29，v4为0.19。在最初的几次迭代中，各节点的权重数值会随着计算发生上下波动。

当执行到第5次与第6次迭代时，数值变化如下 ：

  • 第5次迭代：v1为0.39，v2为0.13，v3为0.28，v4为0.19。
  • 第6次迭代：v1为0.38，v2为0.13，v3为0.28，v4为0.19。可以观察到，此时图中各顶点的权重数值波动已经较小，趋于稳定。当程序监测到相邻两次迭代的数值差小于预设阈值时，即可判定达到收敛并终止循环计算。

基于矩阵的实现

图的拓扑结构可以直接映射为二维的邻接矩阵，从而将复杂的遍历过程转化为线性代数中的矩阵与向量乘法运算。对于包含四个顶点的图结构，其基础邻接矩阵中的元素取值仅包含整数表示边的有无。将基础邻接矩阵进一步转换为概率转移矩阵时，需要根据各源节点的出度，将矩阵列向量中的非零元素进行等比例缩小。

当该转移矩阵与包含各节点当前状态的列向量相乘时，产生的新列向量即为全体节点历经一轮权重传递后的最新状态（做列归一化就相当于上面的除一个∣S(v)∣|S(v)|）：

[tA→AtB→AtC→AtD→AtA→BtB→BtC→BtD→BtA→CtB→CtC→CtD→CtA→DtB→DtC→DtD→D]×[PrAPrBPrCPrD]\begin{bmatrix}t_{A \to A}&t_{B \to A}&t_{C \to A}&t_{D \to A}\\t_{A \to B}&t_{B \to B}&t_{C \to B}&t_{D \to B}\\t_{A \to C}&t_{B \to C}&t_{C \to C}&t_{D \to C}\\t_{A \to D}&t_{B \to D}&t_{C \to D}&t_{D \to D}\end{bmatrix} \times \begin{bmatrix}Pr_A\\Pr_B\\Pr_C\\Pr_D\end{bmatrix}

为了在底层中高效、安全地执行此类张量运算，定义 pagerank 及其需要的 Matrix 类:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
// 头文件保护宏：防止 matrix.h 被多次 #include 导致重复定义
#ifndef PAGERANK_INCLUDE_MATRIX_H
#define PAGERANK_INCLUDE_MATRIX_H

#include <stdio.h>   // printf，用于 print() 输出矩阵
#include <stdlib.h>  // calloc / free，用于连续内存块的分配与释放
#include <assert.h>  // assert，在调试阶段校验前置条件（维度、下标合法性等）
#include <string.h>  // memcpy，用于从外部数组批量拷贝初始数据

// 所有 PageRank 相关数据结构均置于 pagerank 命名空间，避免与标准库或其他模块符号冲突
namespace pagerank {

    /**
     * Row<T> —— 矩阵的「行视图」
     *
     * Matrix 内部并不为每一行单独 malloc，而是持有一块 m×n 的连续缓冲区，
     * 每个 Row 对象仅保存指向该行首元素的指针及列数 n。
     * 这样 matrix[i][j] 的访问路径为：取第 i 个 Row → 取第 j 列元素。
     */
    template <typename T>
    struct Row {
        int n;       // 当前行的列数（与 Matrix 的 n 一致）
        T *rdata;    // 指向该行第一个元素的指针；实际内存由 Matrix::initialize 统一分配
        
        // 默认构造：构造「空行」，尚未绑定任何有效缓冲区
        Row() {
            n = 0;
            rdata = NULL;
        }
        
        /**
         * set_row —— 将 Row 绑定到外部已分配的一段连续内存
         * @param n     列数，必须 > 0
         * @param rdata 指向该行首元素的指针（通常为 buf + i * n）
         *
         * Row 本身不负责释放 rdata；生命周期由 Matrix 管理。
         */
        void set_row(int n, T *rdata) {
            assert(n > 0);
            this->n = n;
            this->rdata = rdata;
        }
        
        /**
         * operator[] —— 列下标访问（语法上写作 row(j) 或 row[j]，取决于运算符重载写法）
         * 返回第 j 列元素的引用，可读写。
         */
        T& operator[](int j) {
            assert(j >= 0 && j < n);
            return rdata[j];
        }
        
        // 析构：仅断开指针引用，不 free(rdata)，避免与 Matrix 双重释放
        ~Row() {
            rdata = NULL;
        }
    };
    
    /**
     * Matrix<T> —— 泛型二维矩阵
     *
     * 内存布局（行主序 row-major）：
     *   buf[0] ... buf[n-1]           ← 第 0 行
     *   buf[n] ... buf[2n-1]          ← 第 1 行
     *   ...
     *   data[i] 指向 buf + i * n
     *
     * 在 PageRank 中，T 通常为 double，矩阵表示转移概率，向量表示各节点当前 rank。
     */
    template <typename T>
    class Matrix {
    protected:
        int m;          // 行数（rows）
        int n;          // 列数（cols）
        Row<T> *data;   // 长度为 m 的 Row 数组，每个元素对应矩阵的一行视图
        
    public:
        // 默认构造：空矩阵，尚未分配任何存储
        Matrix() {
            m = n = 0;
            data = NULL;
        }
        
        /**
         * 参数化构造
         * @param m     行数
         * @param n     列数
         * @param _data 可选，指向 m×n 外部数组；若为非 NULL 则 memcpy 拷贝初始值
         */
        Matrix(int m, int n, const T *_data = NULL) {
            initialize(m, n, _data);
        }
        
        /**
         * initialize —— 核心分配逻辑：一次性申请 m×n 元素缓冲区 + m 个 Row 视图
         *
         * 步骤：
         *   1. calloc 分配 m*n*sizeof(T) 字节，元素默认零初始化
         *   2. 若传入 _data，按行主序 memcpy 整块拷贝
         *   3. new Row<T>[m]，第 i 个 Row 绑定到 buf + i*n
         *
         * 注意：重复调用前应先 clear()，否则旧 buf 会泄漏。
         */
        void initialize(int m, int n, const T *_data) {
            assert(m > 0 && n > 0);
            this->m = m;
            this->n = n;
            T *buf;
            // calloc 相比 malloc 会将内存清零，适合新建矩阵默认全 0
            assert((buf = (T*)calloc(m*n, sizeof(T)))!= NULL);
            if (_data) {
                // 外部 _data 须按行主序排列：elem(i,j) = _data[i*n + j]
                memcpy(buf, _data, sizeof(T) * m * n);
            }
            assert((data = new Row<T>[m])!= NULL);
            // 为每一行建立视图：第 i 行首地址 = 整块 buf 偏移 i*n 个元素
            for(int i = 0; i < m; ++i) {
                data[i].set_row(n, buf + i * n);
            }
        }
        
        /**
         * 拷贝构造函数 —— 深拷贝 other 的全部元素
         *
         * 1. 读取 other 维度；若为空矩阵则直接返回
         * 2. initialize 分配新缓冲区（不拷贝 other 的指针）
         * 3. 双重循环逐元素赋值 data[i][j] = other[i][j]
         */
        Matrix(const Matrix<T> &other) {
            if (&other == this) return;
            this->m = other.row_size();
            this->n = other.col_size();
            if(m <= 0 || n <= 0) return;
            initialize(m, n, NULL);
            for(int i = 0; i < m; ++i) {
                for(int j = 0; j < n; ++j) {
                    data[i][j] = other[i][j];
                }
            }
        }
        
        /**
         * clear —— 释放矩阵占用的堆内存
         *
         * 所有 Row 共享同一块 buf（挂在 data[0].rdata），只需 free 一次；
         * 随后 delete[] data 释放 Row 对象数组。
         * 调用后 data 置 NULL，避免悬空指针。
         */
        void clear() {
            // 所有 Row 共享 data[0].rdata 指向的同一块 buf，只需 free 一次
            if (data != NULL && m > 0 && data[0].rdata != NULL) {
                free(data[0].rdata);
            }
            if (data) {
                delete[] data;
            }
            data = NULL;
            m = n = 0;
        }
        
        // 析构函数：委托 clear 完成资源回收（RAII）
        ~Matrix() {
            clear();
        }
        
        /**
         * operator[] —— 行下标访问，返回第 i 行的 Row 引用
         * 配合 Row 的列访问，实现 matrix[i][j] 二维语法。
         */
        Row<T>& operator[](int i) const {
            assert(i >= 0 && i < m);
            return data[i];
        }
        
        int row_size() const { return m; }  // 返回行数
        int col_size() const { return n; }  // 返回列数
        
        /**
         * operator* —— 矩阵乘法（经典 O(m·n·k) 三重循环）
         *
         * 前置条件：this->n == other.row_size()（左矩阵列数 = 右矩阵行数）
         * 结果维度：m × other.col_size()
         *
         * 计算 res[i][j] = Σ_k ( A[i][k] * B[k][j] )
         *
         * PageRank 场景示例：
         *   转移矩阵 T (n×n) × rank 列向量 v (n×1) → 新一轮 rank (n×1)
         *   即 v' = T * v，对应「每个节点从入链邻居按转移概率汇聚权重」。
         */
        Matrix operator* (const Matrix<T> &other) {
            assert(n == other.row_size());
            int res_col = other.col_size();
            Matrix<T> res(m, res_col);  // 结果矩阵，calloc 已零初始化
            for(int i = 0; i < m; ++i) {
                for(int j = 0; j < res_col; ++j) {
                    res[i][j] = 0;
                    // 内积：第 i 行与 other 的第 j 列做点积
                    for(int k = 0; k < n; ++k) {
                        res[i][j] += data[i][k] * other[k][j];
                    }
                }
            }
            return res;  // 按值返回，可能触发拷贝（取决于编译器 RVO/NRVO）
        }
        
        /**
         * operator= —— 赋值运算符（深拷贝语义）
         *
         * 1. 自赋值检测：&other == this 则直接返回
         * 2. 若当前矩阵已有数据 (m!=0)，先 clear 释放旧内存
         * 3. 按 other 维度重新 initialize，再逐元素拷贝
         *
         * 实现 copy-and-swap 的简化版：先释放再重建，保证维度可变。
         */
        void operator= (const Matrix<T> &other) {
            if (&other == this) return;
            if (this->m!= 0) { clear(); }
            this->m = other.row_size();
            this->n = other.col_size();
            if(m <= 0 || n <= 0) return;
            initialize(m, n, NULL);
            for(int i = 0; i < m; ++i) {
                for(int j = 0; j < n; ++j) {
                    data[i][j] = other[i][j];
                }
            }
        }
        
        /**
         * operator+= —— 同型矩阵逐元素累加（原地修改）
         *
         * 要求两矩阵行数、列数完全一致；否则静默 return（不抛异常）。
         * PageRank 迭代中可用于：rank += damping项、或累加多个转移贡献。
         *
         * 执行 data[i][j] += other[i][j]，即 C += B 的 element-wise 加法。
         */
        void operator+= (const Matrix<T> &other) {
            if (m!= other.row_size() || n!= other.col_size()) return;
            assert(m > 0 && n > 0 && data!= NULL);
            for(int i = 0; i < m; ++i) {
                for(int j = 0; j < n; ++j) {
                    data[i][j] += other[i][j];
                }
            }
        }
        
        /**
         * print —— 按行打印矩阵，元素以 %f 格式输出（假定 T 为浮点类型）
         * 调试用途：可视化转移矩阵或 rank 向量。
         */
        void print() {
            for(int i = 0; i < m; ++i) {
                for(int j = 0; j < n; ++j) {
                    printf("%f ", data[i][j]);
                }
                printf("\n");
            }
        }
    };
}
#endif

基于实现好的 Matrix，即可实现各个网页的 PageRank 更新。

矩阵实现的并行方法

PageRank 每一轮迭代本质是 R' = H × R（H 为转移矩阵，R 为 rank 列向量）。矩阵乘法的各行相互独立，因此可按行划分给多个线程并行计算，最后合并到结果矩阵 res。并行方法的实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#include <utility>   // std::pair, std::make_pair
#include <thread>    // std::thread，创建工作线程
#include <vector>    // std::vector，保存 thread 对象以便 join
#include <cassert>   // assert，校验分区参数合法性
#include <cstdio>    // printf / fprintf
#include "matrix.h"  // 上文定义的 pagerank::Matrix
using namespace pagerank;

// 迭代轮数：重复执行 R ← H×R 共 iters 次
static const int iters = 5;

/**
 * H_arr —— 6×6 转移矩阵的扁平化存储（行主序）
 * 元素 H_arr[i * col1 + j] 表示从节点 j 转移到节点 i 的概率（列随机）
 */
static const float H_arr[] = {
    0, 0, 0.25, 0, 0, 1,
    1, 0, 0.25, 0, 0, 0.2,
    0, 0, 0, 0.5, 0, 1,
    0, 0, 0.25, 0, 1, 0.2,
    0, 1, 0.25, 0.5, 0, 0.5,
    1, 0, 1, 0, 1, 0
};

/**
 * get_range —— 将 [0, total_len) 均分为 partitions 段，返回第 sub_id 段的起止下标
 *
 * @param sub_id      当前分区编号，范围 [0, partitions)
 * @param partitions  分区总数（通常等于线程数）
 * @param total_len   待切分的总行数（此处为 H 的行数 row1）
 * @return            pair<start, end>，左闭右开区间 [start, end)
 *
 * 最后一区（sub_id == partitions - 1）的 end 强制为 total_len，
 * 从而吸收 total_len 不能被 partitions 整除时的余数行。
 *
 * 示例：total_len=6, partitions=3 → 三个区间 [0,2), [2,4), [4,6)
 */
std::pair<int, int> get_range(int sub_id, int partitions, int total_len) {
    assert(partitions > 0 && total_len > 0);
    assert(sub_id >= 0 && sub_id < partitions);
    int each_siz = total_len / partitions;       // 每区至少分配的行数（整除部分）
    int start = sub_id * each_siz;
    int end = (sub_id + 1) * each_siz;
    end = (sub_id == partitions - 1) ? total_len : end;
    return std::make_pair(start, end);
}

int main(int argc, char* argv[]) {
    // H：6×6 转移矩阵；R：6×1 rank 列向量
    const int row1 = 6, col1 = 6;
    Matrix<float> H(row1, col1, H_arr);
    
    // R_arr 须声明为数组；原写法 `float R_arr = {...}` 是非法语法
    static const float R_arr[] = {0.2, 0.2, 0.2, 0.2, 0.2, 0.2};
    const int row2 = 6, col2 = 1;
    Matrix<float> R(row2, col2, R_arr);
    
    // partitions：并行度，每个分区由一个线程负责 H 的若干行与 R 的乘法
    const int partitions = 3;
    if (partitions > row1) {
        fprintf(stderr, "partition is over rows range\n");
    }
    assert(partitions <= row1);
    
    for (int i = 1; i <= iters; ++i) {
        // res 存放本轮 H×R 的完整结果，各线程写入互不重叠的行区间，无需加锁
        Matrix<float> res(row1, col2);
        std::vector<std::thread> threads;
        threads.reserve(partitions);
        
        for (int th_i = 0; th_i < partitions; ++th_i) {
            /**
             * 线程任务（按值传入 th_i，避免捕获循环变量导致的经典并发 bug）：
             *   1. 用 get_range 确定本线程负责的全局行号 [start, end)
             *   2. 从 H_arr 截取对应行块，构造子矩阵 sub1（rows × col1）
             *   3. sub_res = sub1 × R，得到本分区行的局部乘积结果
             *   4. 将 sub_res 写回 res 的全局行 [start, end)
             *
             * H_arr + col1 * start 指向第 start 行首元素（行主序偏移）；
             * Matrix 构造函数会 memcpy 连续 rows×col1 个元素，语义正确。
             *
             * [&] 按引用捕获 H_arr、R、res、col1、col2、row1 等只读或分区写入变量。
             */
            threads.emplace_back([&](int tid) {
                std::pair<int, int> l_rows = get_range(tid, partitions, row1);
                int rows = l_rows.second - l_rows.first;
                Matrix<float> sub1(rows, col1, H_arr + col1 * l_rows.first);
                Matrix<float> sub_res = sub1 * R;
                for (int r = l_rows.first; r < l_rows.second; ++r) {
                    for (int c = 0; c < col2; ++c) {
                        res[r][c] = sub_res[r - l_rows.first][c];
                    }
                }
            }, th_i);
        }
        
        // 等待所有分区线程完成后再读取 res，保证本轮结果完整
        for (int t = 0; t < partitions; ++t) {
            threads[t].join();
        }
        
        printf("iter %d res is:\n", i);
        res.print();
        R = res;  // 深拷贝：下一轮迭代的输入向量更新为本轮输出
    }
    return 0;
}

基于图的实现

使用邻接矩阵的实现方式很直观，但是维护状态转移矩阵存在一个问题：邻接矩阵在实际情况下是稀疏的，即矩阵内部有大量的 0 值元素，存储这些元素并将它们参与计算会消耗大量的系统资源，此时使用基于图的实现更能够节省计算资源。图实现的伪代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// PR_F 处理一条入边 (u, v)，把源点 u 的 rank 贡献累加到目标点 v。
procedure PR_F(u, v)
	// u 的贡献按出度平均分给所有出边终点。
	v.rank += u.rank / out_deg(u)

// VertexMap 在每轮入边贡献累加后应用阻尼系数。
procedure VertexMap(u)
	// 1 - d 是随机跳转项，d * u.rank 是链接转移项。
	u.rank = 1 - d + d * u.rank

procedure PageRank(G=(V,E,D))  // 主循环：输入为包含边，点及其 PR 的图 G
	while i < iters // 迭代终止次数，同时循环终止判断也会考虑迭代是否稳定
		for each v in V // 遍历所有节点
			for each ngh u that satisfies (u, v) in E // 每一次找到一个入边
				PR_F(u, v) // 提取 u 的权重更新 v 的权重
			VertexMap(v) // 使用阻尼系数修正点的 PR
		i += 1 // 进入下一轮迭代

伪代码有一个问题需要注意，在一轮迭代中，每一个点读取到的其他点的 rank，都是上一轮结束时的 rank。实现时需要维护旧 rank 值和本轮更新产生的 rank 值，从而确保计算新值的过程中不会使用其他点在这轮迭代更新的值。

为了兼顾空间效率与访问性能，利用一维数组组合对二维图拓扑进行压缩表达。其中压缩稀疏列方法针对快速提取目标节点的入边集合进行优化。其构建流程如下：

  1. 系统扫描所有原始边缘记录，将边按照目标节点的内部编号升序排列（该例子里就是字典序）。
  2. 按照排序后的顺序，将源节点的标识符连续存储为一个长度与图边数相等的数组（就是中图右列→\to 源点数组）。
  3. 系统统计每一个目标节点所拥有的入边总数量，形成长度等于顶点数量的入度数组。
  4. 生成描述源点范围的偏移量数组，即源点范围数组。

这样的实现极大压缩了存储图信息占用的空间。当需要读某个节点的 rank 时，就进入入度数组按字典序查找；当需要从某个点出发遍历源点时，就在源点数组中按照偏移量逐个访问。

为了进一步提高访问的效率，可以从去除冗余访问入手。在之前伪代码的每一次迭代中，我们对所有点都做更新，试想：如果某一个点提前进入了收敛状态，那么后续基于其做的更新操作都是多余的。因此可以加一次判断，每次把是否 active 的信息 “Pull” 过来进行判断：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// PR_F 只传播源点 u 尚未收敛的 rank 增量。
procedure PR_F(u, v)
	// delta 表示 u 相比上一轮仍需传播的变化量。
	v.rank += u.rank.delta / out_deg(u)

// Pull 版本 PageRank：目标点 v 主动检查入邻居 u 是否活跃。
procedure PageRank(G=(V,E,D))
	while i < iters
		for each v in V 
			for each ngh u that satisfies (u, v) in E
			    if u is active // 判断是否活跃/收敛
					PR_F(u, v)
			// 汇总所有活跃入邻居贡献后，更新 v 的 rank 和活跃状态。
			VertexMap(v)
		i += 1

此时 PR_F 不是直接用 rank 做加法，v.rank 也不会在新一轮迭代开始时重置成 0，而是基于增量做更新，如果按照上一版伪代码的 PR_F 规则计算是完全错误的！这种方法也有一个缺点——每次都做 active 判断会拖慢算法的运行速度。要优化这个缺点，可以从源节点出发，让每一个顶点的 rank 在发生更新时，主动把变化 “Push” 给指向的节点，从而避免了冗余的判断：

1
2
3
4
5
6
7
8
9
10
11
12
13
// Push 版本 PageRank：活跃源点主动向出邻居传播增量。
procedure PageRank(G = (U, E, D)) 
	while i < iters 
		// 活跃节点主动把 Delta Push 出去
		for each active vertex u in U 
			for each ngh v that satisfies (u, v) in E 
				PR_F(u, v) // 将 u 的 delta 累加到 v 的接收缓冲区 (v.temp_delta)

		// 统一结算本轮收到的所有数据
		for each v in V // 或者仅遍历收到了消息的 v
			VertexMap(v) // 执行 v.rank += v.temp_delta，并判断是否需要变为 active
			// 结算后清空 v.temp_delta，供下一轮接收新的增量。
		i+=1

图实现的并行方法

对于矩阵的并行化很简单——把状态转移矩阵每一行与 rank 矩阵的运算分摊到多个线程上即可。而使用基于图的方法时，并行方法就没有矩阵那么直观。和之前最小生成树的并行一样，在图实现中将图进行划分发给不同的线程进行更新。首先给定一个原始图：

545

为了简化运算，只将图划分给两个线程，分别处理 PointSet1{1,2,3}PointSet_1\{1, 2, 3\} 和 PointSet2{4,5,6}PointSet_2\{4, 5, 6\} 两个点集。对于两个点集，划分出来的子图中保留对应点集中点之间的边以及这些点指向的边（因为采用 active 优化的情况下需要同步更新指向的点）：

划分完成后，对于子图分别做 PageRank 处理。需要注意的是，这种方法可能出现两个线程同时写一个 rank 的情况（例如线程 1 在将顶点 1 的数据 Push 给顶点 2 时，线程 2 同时 将顶点 6 的数据 Push 给了顶点 2），所以具体实现时需要做好同步和加锁处理。

上面的实现基于两个点集的出边进行计算，在这种方法中基于两个点集的入边进行计算：

893

这种情况下，两个子图中的点不可能同时将数据 Push 给同一个点。可以这样考虑：

A∈Graph1A \in Graph_1，分两种情况：

  1. A∈PointSet1A \in PointSet_1，那么 A∉PointSet2A \notin PointSet_2，它在子图 2 中只有出边，那么子图 2 中一定没有边指向 AA；
  2. A∈PointSet2A \in PointSet_2，它在子图 1 中只有出边，那么子图 1 中一定没有边指向 AA。

同理也可以证明在 A∈Graph2A \in Graph_2 的情况下，也不被两个点同时 Push。这样就不需要额外的同步操作来防止数据冲突了。所有顶点的当前数据被连续的存储在内存数组中，线程 1 和线程 2 可能同时读取某个顶点的当前数据，但并不会修改该数据。如果每个线程都按照存储顺序依次访问各个顶点，便可以最大程度的减少 cache miss 次数，从而有效提高性能。