目录#

共享内存的操作#

Windows下可以用Windows的API进行操作，而Linux下则可以通过POSIX的API进行操作，虽然过程有点区别，但终究是大同小异。

创建&数据写入#

为了简单测试，这里就直接用Windows版Python自带的mmap当实验：

1
2
3
4
5
6
7
8
9
10
11
12
13

>>> import numpy as np
>>> a = np.random.randn(32*4096,4096)  # 先占用4G内存
>>> a[0,0]
-0.5271477716515403
>>> a[-1,-1]
1.1194140920933267
>>> b = a.tobytes()  # 额外的4G内存，peak 8G
>>> del a  # 释放4G内存
>>> import mmap
>>> m = mmap.mmap(0, 32*4096*4096*8, 'Local\\MyShMem1')
>>> m.write(b)  # 额外的4G内存，peak 8G
>>> del b  # 释放4G内存
>>>

这段代码会常占4G的内存，顶峰的内存会达到数据量的2倍。（在写入共享内存的时候）

注：Linux版Python不支持这种创建共享内存的方式。

数据读取#

数据读取是用matlab的mex写的：

shmem_win_data_read.c：

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

#include <mex.h>
#include <stdio.h>
#include <windows.h>
#pragma comment(lib, "user32.lib")

TCHAR szName[] = TEXT("Local\\MyShMem1");
void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    //mxUint64 ptr = mxGetUint64s(prhs[0])[0];
    mxArray* pArr = mxCreateDoubleMatrix(32*4096, 4096, 0);
    HANDLE hMapFile = OpenFileMapping(FILE_MAP_ALL_ACCESS, FALSE, szName); // param: R/W, do not inherit the name, name of mapping object
    if (hMapFile == NULL) {
        printf("Failed to open file mapping: %d\n", GetLastError());
        return;
    }
    double* b = (double*)MapViewOfFile(hMapFile, FILE_MAP_ALL_ACCESS, 0, 0, 32*4096*4096*sizeof(double)); // Last: bufsize
    if (b == NULL) {
        printf("Failed to map view of file: %d\n", GetLastError());
        return;
    }
    for(int i = 0; i < 5; i++)
        printf("%lf\n", b[i]);
    double* originalArr = mxGetPr(pArr);
    mxFree(originalArr);
    printf("Original ptr: %lld\n", originalArr);
    mxSetPr(pArr, b);
    printf("Current ptr: %lld\n", mxGetPr(pArr));
    plhs[0] = pArr;
    if (nlhs > 1) {
        // handle
        mwSize d[] = { 1 };
        pArr = mxCreateNumericArray(1, d, mxUINT64_CLASS, mxREAL);
        *(unsigned long long*)mxGetData(pArr) = (unsigned long long)hMapFile;
        plhs[1] = pArr;
    }
}

然后编译，运行：

1
2
3
4
5

mex 'shmem_win_data_read.c' -g
a = cell(1);
[a{1}, h] = shmem_win_data_read;
b = a{1};
% 用变量b进行正常读写

返回的第一个参数是数据，至于为什么要用cell，就是为了规避matlab自带的GC机制。一旦将数组的指针利用mxSetPr设置为一个非通过调用mxMalloc得到的值（如c语言自带的malloc），GC机制尝试对它进行free的时候，后果自己可以想象。

返回的第二个参数是共享内存的handle，需要在使用完后释放掉。

cell的内存机制后面再说，注意的是，在调用mxSetPr之前，一定要用mxGetPr和mxFree把最初创建返回数组pArr时所申请到的内存释放掉，否则就会发生内存泄漏。

可以看出，数据是完全一致的。

资源释放#

一旦用完，就该妥善处理后事了（笑）。你大可以试试clear直接清掉，看看matlab会不会崩。

shmem_win_data_free.c：

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

#include <mex.h>
#include <stdio.h>
#include <windows.h>
#pragma comment(lib, "user32.lib")
#include <matrix.h>
#include <stdlib.h>

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    mxArray* cell_arr = mxGetCell(prhs[0], 0);
    double* new_ptr = (double *)mxMalloc(sizeof(double) * 5);
    double* original_ptr = mxGetPr(cell_arr);
    new_ptr[0] = 1.2;
    new_ptr[1] = 0.0;
    new_ptr[2] = 0.0;
    new_ptr[3] = 0.0;
    printf("new ptr: %lld\n", new_ptr);
    printf("original ptr: %lld\n", original_ptr);
    mxSetPr(cell_arr, new_ptr);
    printf("current ptr: %lld\n", mxGetPr(cell_arr));
    mxSetM(cell_arr, 2);
    mxSetN(cell_arr, 2);
    UnmapViewOfFile(original_ptr);
    if (nrhs > 1) {
        HANDLE hMapFile = *(HANDLE*)mxGetData(prhs[1]);
        CloseHandle(hMapFile);
    }
}

为了简单实验起见，这里就直接用把内存改回matlab的托管内存了（用mxMalloc申请的内存），然后手动把之前改的内存释放掉。

1
2

mex 'shmem_win_data_free.c' -g
shmem_win_data_free(a, h);  % 这里就用a作为参数传入

合理利用matlab的“特性”，a{1}和b都会同时修改为一个新的数组，接着就可以直接clear掉，不需要再担心什么了。

完整使用#

1
2
3
4
5
6
7
8
9

v = zeros(1, 4096);
parfor x = 1:4096
    a = cell(1);
    [a{1}, h] = shmem_win_data_read;
    b = a{1};
    c = sum(b(:, x));
    v(x) = c;
    shmem_win_data_free(a, h);
end

这时候大可安心看系统的内存使用情况，妈妈再也不用担心内存占用啦。

番外：直接使用mxSetPr的效果#

mxSetPr在mex中一般用法效果如下：

1
2
3
4
5
6
7
8
9
10
11
12

#include <mex.h>
#include <stdio.h>

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    const mxArray* arg = prhs[0];
    double* ptr = mxGetPr(arg);
    printf("Original dataArray PTR: 0x%p\n", ptr);
    double* new_ptr = (double*)mxMalloc(sizeof(double) * mxGetM(arg) * mxGetN(arg));
    printf("Allocated dataArray PTR: 0x%p\n", new_ptr);
    mxSetPr(arg, new_ptr);
    printf("New dataArray PTR: 0x%p\n", mxGetPr(arg));
}

上面代码的作用：直接申请一块新的内存，将这块内存的地址直接赋值到matlab（作为输入参数传入mex函数）的矩阵数组中。运行结果如下：

可以看到，在mex函数体内对数组进行的任何修改是无法直接影响到输入参数的。

番外2：必须用cell中的元素作为输出变量#

还有一个关键点，就是必须用a{1}当输出变量，如：

1
2
3
4
5
6
7

a{1} = randn(5);
b = a{1};
% 通过b进行数据读取

do_cleanup_within_cell(a);
% 上面的函数能够同时作用于a（输入变量）和b（通过a引用赋值）
% 执行结果：a = 1×1 cell, a{1} = [], b = []

而把一般变量作为输出，再套用cell，则只能影响到cell里面的元素，不会影响到最先赋值的变量

1
2
3
4
5
6
7
8

t = randn(5); % 先赋值到一个新变量
a{1} = t; % 再通过引用赋值给cell
b = a{1};
% 通过t或者b进行数据读取

do_cleanup_within_cell(a);
% 上面的函数仅能作用于a（输入变量）和b（通过a引用赋值），而最原始的t则不受影响
% 执行结果：t = 5×5 double, a = 1×1 cell, a{1} = [], b = []

其中do_cleanup_within_cell可为如下（直接把输入数组清空为[]）：

1
2
3
4
5
6
7
8

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    mxArray* cell_arr = mxGetCell(prhs[0], 0);
    double* original_ptr = mxGetPr(cell_arr);
    printf("original ptr: 0x%p\n", original_ptr);
    mxSetM(cell_arr, 0);
    mxSetN(cell_arr, 0);
    mxSetPr(cell_arr, NULL);
}

反正个人猜测这种现象肯定跟matlab的GC是脱不了钩的，由于没有官方的内存介绍文档，所以也不深究了。写代码讲究一件事，那就是能用就行。

数据类型#

除了double、float以及(u)int8/16/32/64外，还有complex、sparse、cell、struct和object等数据类型，一个完整的check如下：

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

#include <mex.h>
#include <matrix.h>
#include <stdio.h>

#define CHECK_FUNC(x) printf(#x": %d\n", x(prhs[0]))

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    // 1. mxArray attributes
    CHECK_FUNC(mxIsNumeric);
    CHECK_FUNC(mxIsComplex);
    // 2. create, query, and access data types
    // 2.2 noncomplex float
    CHECK_FUNC(mxIsScalar);
    CHECK_FUNC(mxIsDouble);
    CHECK_FUNC(mxIsSingle);
    // 2.3 noncomplex integer
    CHECK_FUNC(mxIsInt8);
    CHECK_FUNC(mxIsUint8);
    CHECK_FUNC(mxIsInt16);
    CHECK_FUNC(mxIsUint16);
    CHECK_FUNC(mxIsInt32);
    CHECK_FUNC(mxIsUint32);
    CHECK_FUNC(mxIsInt64);
    CHECK_FUNC(mxIsUint64);
    // 2.6 sparse
    CHECK_FUNC(mxIsSparse);
    // 2.8 character
    CHECK_FUNC(mxIsChar);
    // 2.9 logical
    CHECK_FUNC(mxIsLogical);
    // 2.10 object
    // mxIsClass(mxArray* arr, char* className)
    // 2.11 structure
    CHECK_FUNC(mxIsStruct);
    // 2.12 cell
    CHECK_FUNC(mxIsCell);
}

接下来简单地验证完Linux的共享内存之后，就可以开始造（已经重复）的轮子了。

Linux版的碎碎念#

手头上只有Win版的R2018b和Linux版的R2017b服务器，想着版本差距不算大，直接改改代码，把WinAPI的Named shared memory改成POSIX的shm_open大概就行了。

把代码扔过去，编译，运行，哦豁完蛋，报了个在Windows下未曾出现过的崩溃错误：

突然告诉我还有个16Byte的header？喂不是吧大哥，我只把数据扔共享内存上了你告诉我header？

好吧，既然你想读，那就让你读呗。我把共享内存创建大一点，把数据扔后面一点，前面留16Byte的零，总算可以吧。然后：

你塔喵还要查header是吧？好好，我把mxGetPr往前16个字节的数据一并复制过去总行了吧。

嗯，这次好像真行了，而且能够正常退出没崩。

至于header是什么，说实话自己也没弄懂。随便弄了几个矩阵，往前挪16个字节print一下，结果如下：

131072*4096大小的double矩阵的header：
00 00 00 00 01 00 00 00 CE FA ED FE 20 00 10 00

7*65537大小的double矩阵：
40 00 38 00 00 00 00 00 CE FA ED FE 20 00 20 00

7*65538大小的double矩阵：
70 00 38 00 00 00 00 00 CE FA ED FE 20 00 10 00

7*65538大小的int32/uint32/single矩阵：
40 00 1C 00 00 00 00 00 CE FA ED FE 20 00 10 00

前面8字节明显跟数组的大小相关，那么后面8字节呢，跟type也不相关啊……特别是对sparse、complex这种数据类型也不变，而且倒数第二个0x10偶尔会变成0x20，也毫无规律可言。

目前只能推测跟matlab的内部GC有关（因为改一下清空变量的时候就会崩掉）。

（挖坑，目前只是试验可行性，不一定100%能做出什么东西出来）

目录#

引言#

无论跑什么代码也好，只要数据的访问是相互独立，没有任何依赖关系的，一般情况下都可以引入并行计算增加CPU/GPU的利用率。matlab也不例外，最直接傻瓜的：将for改成parfor，就可以利用matlab自带的并行计算工具提高代码的运行效率。

举个栗子，就用手头上推荐系统的NDCG计算，每个用户的NDCG是相互独立的，所以为每个用户计算NDCG这一个过程就可以通过改成parfor i = 1:n实现并行化，最后再将结果取平均。

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

function ndcg = metric_ndcg(R, R_hat, k)
n = size(R, 1);
ndcg = zeros(1, n);
invalid_users = zeros(1, n);
RT = R';
RT_hat = R_hat';
parfor i = 1:n
    Ri = RT(:, i); % [1, items]
    Ri_hat = RT_hat(:, i)';
    non_empty_entries = find(Ri ~= 0);
    Ri = Ri(non_empty_entries);
    Ri_hat = Ri_hat(non_empty_entries);
    [~, idx] = sort(Ri_hat, 'descend');
    dcg_Ri = Ri(idx);
    dcg = dcg_at_k(dcg_Ri, k, true);
    
    idcg_Ri = sort(Ri, 'descend');
    idcg = dcg_at_k(idcg_Ri, k, true);
    if idcg ~= 0
        ndcg(i) = dcg / idcg;
    else
        invalid_users(i) = 1;
    end % if
end % for i
invalid_users = sum(invalid_users);
if n ~= invalid_users
    ndcg = sum(ndcg, 2)' / (n - invalid_users);
else
    ndcg = nan;
end
end % function

function dcg = dcg_at_k(r, k, p)
n = min([k length(r)]);
d = log2(2:n+1)';
rn = r(1:n);
if p
    rn = power(2, rn) - 1;
end
dcg = sum(rn ./ d);
end % function

当然，跑小的数据集不成问题，只不过评分矩阵R一旦大起来，内存的使用就会直线飙升。为什么？多亏了matlab的沙雕并行计算的内存机制，下面一一细数它的罪恶。

matlab的内存机制#

自己接触matlab不算多，平时也就写写代码，能跑就行的那种。但是，偶尔为了优化一下代码，摸点底层机制还是必要的。

为了探究这个机制，需要用到matlab中的mex，它是matlab提供的能调用其他编程语言的一个工具。这里就用c语言作为栗子，毕竟有它就能够直接实现很多骚操作：

官方文档：https://www.mathworks.com/help/matlab/matlab_prog/memory-allocation.html

mex c：从入门到治疗高血压#

首先需要#include <mex.h>，matlab调用c语言的函数void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])也是不可少的，其重要性相当于平时的main()。nlhs和nrhs两个值代表输入和输出的参数个数，输入的参数都保存在prhs数组中，而输出的参数则需要自己创建，赋值到plhs数组中。

这个栗子只是简单输出一下第一个输入参数的内存位置，函数体内第一行的mxArray PTR代表目前作为参数的这个矩阵对象的内存地址，而第二行dataArray PTR则是用mxGetPr获取这个对象里面的实际数据的内存位置，第三行就简单输出第一个数据的值。

编译和调用过程如下：

这两个指针的含义可以用一段等效的c++代码直观反映出来，在这个栗子中，mxArray PTR相当于&a，而dataArray PTR相当于a.array：

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

#include <string.h>

class mxArray {
public:
    double* array; // 实际保存数据的指针
    mxArray (int size_1, int size_2) { // 构造函数
        // ...
    }
    ~mxArray () {
        // ...
    }
};

// 调用的mxGetPr的等效代码
double* mxGetPr(const mxArray* arr) {
    return arr->array;
}

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    // ...
}

int main() {
    mxArray a(32*4096, 4096); // a = zeros(32*4096, 4096)
    const mxArray* b[] = { &a }; mexFunction(0, NULL, 1, b); // double_ptr(a)
    return 0;
}

矩阵的存储方式#

与python中的numpy不同，matlab是按维度索引从左到右进行存储的，比如a=randn(3,3)，在内存上则以a(1,1) a(2,1) a(3,1) a(1,2) a(2,2) a(3,2) a(1,3) a(2,3) a(3,3)的顺序分布，所以提高CPU缓存命中率的读取方式是a(:,1)这类固定后面若干维度的范围读取形式。

矩阵的GC机制#

matlab的GC是以数组的引用数为依据的，分配内存地址的时候以及进行读写操作时，matlab会自动跟踪并调整引用数，当引用数为0（意味着已经没有变量用到这段数据）的时候，GC会自动释放掉这段内存。

常见的赋值方式#

首先声明一下，赋值方式那名称都是我瞎起的……

为了更能掌握其中的数组指针变化，我写了一个新的mex函数用来创建数组：

第一种是最常见的直接赋值，即将右边的数组直接赋值到变量上：

1
2

a = new_arr; % <--
a = new_arr; % <--

并且第二条赋值语句完成后，变量a的数组指针地址会发生变化。

这种赋值操作的具体顺序是（当然不是每步都有验证，也没确定对不对）：

1. zeros创建一个数组，matlab为其分配新内存，将引用数记为1，全部置0后返回这个数组
2. 创建变量a，将a中的数据指针的地址赋值为上面的数组地址，增加引用数（现在是2）
3. 清理调用zeros函数时创建的所有数组，包括减少zeros创建并返回的数组的一个引用数（现在是1），这条清理规则如果熟悉c++对象生命周期的话应该不难理解
4. randn创建一个数组，matlab为其分配新内存，将引用数记为1，赋值完成后返回这个数组
5. a变量的数据指针重新指向新的数组地址，新的数组的引用数增加为2，同时减少原本旧的数组的引用数（变为0）
6. 清理调用randn函数时创建的所有数组，这时将新数组的引用数减少为1
7. GC检测到有引用数为0的数组，释放原来那个数组的内存

第二种也是很常见的引用赋值，即将右边的数组没做任何修改直接赋值到新的变量上：

1
2

a = randn(32*4096, 4096);
b = a; % <--

赋值完成后，变量b的数组指针与变量a一致。（同时这段数组的引用数会自增）至于说数据修改之后怎么影响，这就是matlab的另外一个机制：写入复制（copy on write）的事了。

第三种，一般很少用，我叫它覆盖赋值，前提是等号两边的数组个数要一致。赋值后，右边的数组数据会复制到左边数组中，不会改变左边的数组指针的地址：

1
2

a = new_arr;
a(:) = new_arr; % <--

写入复制机制#

当数组的引用数大于1时，只要对其中一个引用了该数组的变量进行写入，都会触发写入复制机制，另外创建一个数组保存这个修改，原数组引用数-1。而引用数只有1的时候，则会写入到原数组地址中。

栗子：

1
2
3
4

a = new_arr;
b = a; % 引用赋值
a(1) = 233; % <--
b(1) = 666; % <--

然而有一个很有趣的“特性”：mex不受这个机制的约束，也就是说，只要拿到地址，就可以绕过这个机制进行写入，如：

当然，还有更绝的：

matlab并行计算的内存机制#

举个栗子，假如用E5的CPU跑parfor，数据矩阵A占4G，matlab会创建12个并行的计算进程占满12个CPU核心。整个架构管它叫master-slave也好，server-client也罢，反正master/server会向slave/client发送进行计算所需要的所有数据，slave/client独立进行计算，并将结果返回到master/server中。

那么沙雕matlab会怎么操作呢？

把数据矩阵A再复制12遍。没错，是再复制12遍。一个4G，那么加起来13个就会占掉52G的内存。如果看到这里已经感到不适，该去医院检查高血压了。

首先来几个FAQ：

Q1：既然matlab能调用c/c++代码，那我在c/c++代码里面手动创建线程，再调用matlab，不用它parfor不就行了吗？
A1：No，当你在c/c++创建新线程再调用matlab的时候，它已经注定要crash了。

Q2：那我在每次循环中只处理A中的某一行/某一列不就行了吗？
A2：No，就算每次只处理一行/列也好，两行/列也罢，matlab终究会把整个数据矩阵A完完整整地复制一遍的。

Q3：So easy，那我自己造个轮子，将共享内存的地址赋值到matlab，再手动释放，避开matlab的主动复制不就行了吗？
A3：可以，但愿你能从每次操作内存引发的crash中吸取教训。

因此，相对于从matlab自带的进程间通讯机制下手而言，选择相对容易解决的内存机制下手是比较明智的，好好利用matlab的一些“特性”应该不难实现。

To be continued.

需要sudo权限，matlab的安装iso文件，可选的有激活用的.lic许可证文件和一些破解验证的.so文件。

1. 挂载iso：sudo mkdir matlab_iso && sudo mount -o loop R2017b_glnxa64.iso matlab_iso/
2. 复制iso里面的文件：mkdir matlab && cp -R matlab_iso/ matlab/
3. 取消iso挂载：sudo umount matlab_iso/
4. 把配置文件加写入权限（官方建议是复制一份来改的）：cd matlab/matlab_iso && chmod +w . activate.ini installer_input.txt
5. 修改installer_input.txt，主要修改以下几个字段：
   destinationFolder=/usr/local/R2017b指定安装的目录
   fileInstallationKey=xxxxx-xxxxx-....指定激活的那串密钥
   agreeToLicense=yes同意license
   outputFile=/tmp/mathworks.log指定安装log保存的位置
   mode=silent
activationPropertiesFile=activate.ini其实这个改不改应该不影响，毕竟后面要手动激活
6. 安装matlab：sudo ./install -inputFile installer_input.txt（建议用绝对路径避免出现其他麻烦）
7. （可选）覆盖原来的libmwservices.so：cp .../libmwservices.so /usr/local/R2017b/bin/glnxa64/libmwservices.so
8. 修改激活配置文件activate.ini：
   activateCommand=activateOffline指定离线激活
   licenseFile=.../license_standalone.lic指定许可证文件的位置（建议用绝对路径）
9. 激活matlab：/usr/local/R2017b/bin/activate_matlab.sh -propertiesFile activate.ini

大概就差不多了，matlab就在/usr/local/R2017b/bin/matlab里面，往.bashrc结尾加上export PATH=$PATH:/usr/local/R2017b/bin即可。

正好这时间快开无限池终本了，本来还想测试一下刚实现的助战技能等级识别，结果一上来就死在礼装检测上了。无限池那张加成礼装识别不出来。

俺寻思，俺是用SIFT来匹配的，怎么也不应该呀。看看从屏幕截取到的的助战礼装栏：

没有看出任何问题，再看看礼装的完整卡面：

拍拍脑袋都知道，这肯定能匹配上的。为了避免自己在奇怪的地方出错，特意手动跑了一遍SIFT：

依然没有任何问题……那时候真就在怀疑人生。后面拍脑袋才发现，好像图片不对！我爬的是礼装图鉴的缩略图，也就是：

好家伙，这图片能匹配得上就有鬼了。你月也是截图鬼才，默默改代码去了。

附：缩略图通过SIFT与原图匹配的情况：

记录配置OpenWrt的IPv6 NAT和DHCP踩过的一些坑。

Read More