在上一篇文章中，采用了写C模块的方式来加速mandelbrot分形图的计算，虽然比起用python写要快得多，但在运行时观测各个cpu的占用率是就会发现，在没有运行其他高负载任务的情况下，只有1个核的占用率会达到100%，会出现"一核有难，X核围观"的情况。尤其是现代的主流计算机，基本上都是多核心的，计算中只占用一个核，效率是很低的，极大地浪费时间与资源。
使python进行多核并行的方式有很多(原生的multiprocessing、原生的MPI接口mpi4py等)，其中mpi4py方式符合MPI标准，可以方便地在大型多节点分布式集群环境下进行多核并行计算。
MPI(Message Passing Interface)是一套标准，一个跨语言的通讯协议，用于编写并行计算机。在不同的硬件环境与不同的编程语言下有不同的实现，在python中，可以使用mpi4py来实现多核并行计算。
在MPI的早期发展，计算机主要是单核架构，每个核与内存是分布在不同的物理机器上的，所以主要发展的是在不同进程中消息传递。随着现代计算机技术的发展，多核架构的机器成为主流，通常在一个节点上有若干个核共用一块物理内存，相比在进程间消息传递，在同一个物理机上采用共同读写一块内存的方式，可以节省通信所造成的开销，从而提高计算效率。MPI-3版本新增加了一种内存共享机制，进程可以将自己的一部分内存通过窗口共享给其他的进程，共享的内存可以像普通的内存一样直接存取，而不需要通过通信，从而减少内存数据复制所带来的性能损失与由此所占用的多余的内存。
一般来说，在并行模型中可以采用多进程与多线程的，线程之间的内存与变量是共享的，进程之间的内存与变量是独立的，然而对于Python(指Cpython)来说，由于GIL(Global Interpreter Lock，全局解释器锁)的存在，所以在多线程中，在同一时间实际上只会有一个线程在运行。所以多线程并不能提高计算的效率，而多进程由于每一个进程都有自己的内存空间，有自己的GIL，所以可以同时运行。至于为什么会有GIL这个东西，这大概是历史遗留问题，毕竟在刚有python的时候，还没有线程这个概念……(实际上，python和MPI都是1990年左右诞生的)
mpi4py的安装并不复杂，通常情况下用pip安装即可：
pip install mpi4py然后可以测试一下能否正常运行：
mpiexec -n 5 python -m mpi4py.bench helloworld在MPI中，会启动指定数量的进程来同时运行，进程之间的通信与协调则是由mpi4py提供的API来完成的。在python中导入mpi4py的模块mpi4py.MPI,并获取进程总数与进程的序号(rank)：
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
print("Hello World from process %d of %d".format(rank, size))回到例子，画一个mandelbrot图，需要知道在各个点的具体迭代次数，而且各个点之间的计算又是独立的，可以在MPI框架下分发给不同进程，然后每个进程计算自己的部分，最后再合并计算结果。合并计算结果这一步可以通过共享内存来完成，先划出一块共享内存(numpy 数组)，然后每个进程都会把自己的计算结果直接写入共享内存(numpy 数组)，最后直接把这个结果数组用matplotlib画出来即可。
先写划出共享内存，这里可以让0号进程来分配一块连续的内存空间，把这个内存空间的地址分发给其他进程，然后每个进程都会把自己的计算结果写入这块内存空间，所有的进程都依据这一块连续的内存空间生成一个numpy数组，然后把计算结果写到这个数组中，最后再用matplotlib中的imsave保存成png文件即可。
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
# 获取基本信息
itemsize=MPI.SHORT.Get_size()
# 由于迭代次数是一个并不会很大的整数，这里跨域用short来存储，对应的numpy数据类型是int16
X=1920
Y=1080
# 生成图片的尺寸，这里按照1080P的分辨率来设置
if rank == 0:
bits=X*Y*itemsize #0号进程需要分配的内存空间大小
else:
bits=0 # 其余进程不需要分配内存
win=MPI.Win.Allocate_shared(bits, itemsize, comm=comm) # 分配内存空间
buf,isize=win.Shared_query(0) # 获取内存空间(0号进程拿到的)的地址
buf = np.array(buf, dtype='B', copy=False)
ary = np.ndarray(buffer=buf, dtype='int16', shape=(,X*Y)) # 将内存空间转换成numpy数组
这样就拿到了一个在全部线程里都能一起读写的numpy数组了，然后再把整个1920X1080个像素点分发给这些进程计算，然后再把结果写入这个数组中，最后再用matplotlib中的imsave保存成png文件即可。
dx=(x[1]-x[0])/X
dy=(y[1]-y[0])/Y
# x,y为在复平面上坐标的范围，原则上应该是按照16：9的比例来设置的
for item in range(rank,X*Y,size): # 每个进程计算一部分
yy=item//X
xx=item%X
ary[item]=mcssdog(x[0]+xx*dx,y[0]+yy*dy,255,4)# 计算每个像素点的迭代次数,这里的255是迭代次数的上限，这个函数具体内容与作用可以看上一篇文章
comm.Barrier() # 强制同步 即所有进程都运行到这里之后才能继续执行
if rank==0:
plt.imsave('out.png',ary.reshape(Y,X),vmax=255,vmin=0,cmap='jet') # 保存成png文件这样就可以画出一张图来，但是如果需要若干张图用来制作gif动图或者视频，就可以在保存图片的时候也并行处理了，一次性计算完多张图片的数据，然后把保存图像的任务分发给各个进程来完成，这样就可以提高效率了。
下面是完整的程序，可以从一个写有中心点坐标与宽度的txt文件生成出一系列的图片
C 代码mc.c，具体内容可以看上一篇文章 ：
#include <Python.h>
static short mpdog(double a, double b, int T, float M)
{
double x=0;
double y=0;
double c=0;
double d=0;
int i = 0;
for (; i < T; i++)
{
c=x*x-y*y;
d=2*x*y;
c=c+a;
d=d+b;
if (c*c+d*d>M)
{
return i;
break;
} else
x=c;
y=d;
}
return i;
}
static PyObject * mcssdog(PyObject * self, PyObject * args)
{
double a = 0;
double b = 0;
int c=0;
float d=0;
if (!PyArg_ParseTuple(args, "ddif", &a, &b, &c, &d))
return NULL;
short r = mpdog(a, b,c,d);
return Py_BuildValue("h", r);
}
static PyMethodDef accuMethods[] =
{
{"mcssdog", mcssdog, METH_VARARGS, "add two floats"},
{NULL, NULL, 0, NULL}
};
static struct PyModuleDef caadogm =
{
PyModuleDef_HEAD_INIT,
"mcssdog",
"caadogddd",
-1,
accuMethods
};
PyMODINIT_FUNC PyInit_mcssdog(void) {
return PyModule_Create(&caadogm);
}
编译成动态链接库
gcc -I/usr/include/python3.8 -o mcssdog.so mc.c -fPIC -sharedpython部分mat.py
import numpy as np
from mpi4py import MPI
import time
from matplotlib import pyplot as plt
from mcssdog import mcssdog
from genlist import genlist,readpath
from sys import argv
if len(argv)<=1:
filepath="path.txt"
else:
filepath=argv[1]
center_list,r=readpath(filepath)
rant=16/9
X=1920
Y=1080
mMAX=255
FPS=120
Xlist,Ylist=genlist(center_list,r,rant,N=FPS)
Num=len(Xlist)
N=0
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
itemsize=MPI.SHORT.Get_size()
if rank == 0:
bits=X*Y*itemsize*size
else:
bits=0
win=MPI.Win.Allocate_shared(bits, itemsize, comm=comm)
buf,isize=win.Shared_query(0)
buf = np.array(buf, dtype='B', copy=False)
ary = np.ndarray(buffer=buf, dtype='int16', shape=(size,X*Y))
def mand(x,y,X,Y,M):
dx=(x[1]-x[0])/X
dy=(y[1]-y[0])/Y
for item in range(rank,X*Y,size):
yy=item//X
xx=item%X
M[item]=mcssdog(x[0]+xx*dx,y[0]+yy*dy,mMAX,4)
comm.Barrier()
if rank==0:
start_time=time.time()
while N<Num:
huan=N%size
mand(Xlist[N],Ylist[N],X,Y,ary[huan])
if huan==size-1:
comm.Barrier()
filename="pic/mandelbrot{}.png".format(N-(size-1-rank))
plt.imsave(filename,ary[rank].reshape(Y,X),vmax=mMAX,vmin=0,cmap='jet')
N=N+1
else:
S=Num%size
if rank<S:
filename="pic/mandelbrot{}.png".format(N-(S-1-rank))
plt.imsave(filename,ary[rank].reshape(Y,X),vmax=mMAX,vmin=0,cmap='jet')
if rank==0:
end_time=time.time()
print("Time:{}s".format(end_time-start_time))其中genlist.py:
import numpy as np
def readpath(path):
center_list=[]
r=[]
f = open(path, 'r')
for ch in f.readlines():
if ch[0] == '#' or ch[0] == '\n':
continue
else:
p = ch.split()
p = [float(i) for i in p]
center_list.append(p[0:2])
r.append(p[2])
f.close()
return center_list,r
def genlist(p,r,rant,N=30):
"""
生成参数列表
"""
if len(p)<=1 or len(r)<=1:
return None
x=len(p)-1
Xlist=np.zeros([x*N,2])
Ylist=np.zeros([x*N,2])
for i in range(x):
Xlist[i*N:(i+1)*N,0]=np.linspace(p[i][0],p[i+1][0],N)-np.linspace(r[i],r[i+1],N)/2*rant
Ylist[i*N:(i+1)*N,0]=np.linspace(p[i][1],p[i+1][1],N)-np.linspace(r[i],r[i+1],N)/2
Xlist[i*N:(i+1)*N,1]=np.linspace(p[i][0],p[i+1][0],N)+np.linspace(r[i],r[i+1],N)/2*rant
Ylist[i*N:(i+1)*N,1]=np.linspace(p[i][1],p[i+1][1],N)+np.linspace(r[i],r[i+1],N)/2
return Xlist,Ylist
路径文件
# 路径文件
-0.6 0 2.8
0.2775 0.008 2.8
0.2775 0.008 0.1
0.2775 0.0085 0.05
0.2774 0.0086 0.005
0.2773 0.0086 0.001
0.2770 0.0086 0.0005
0.2770 0.0086 0.0001
0.2770 0.0086 0.00005
0.2770 0.0086 0.00001
0.2770 0.0086 0.000005
0.277005 0.0086 0.000005
0.277010 0.0086 0.000005
0.277015 0.0086 0.000005
0.277015 0.008605 0.000005
0.277015 0.00861 0.000005
0.277015 0.008615 0.000005
0.277015 0.00862 0.000005运行
mpirun --use-hwthread-cpus -np 16 python mat.py在i7-10700上运行，使用全部16个核，共画出2041张图片，用时14min。
成果
视频如下
参考资料