曼德勃罗集程序
Zn+1=Z2n+C,
曼德勃罗特集是易并行计算
任务分配和动态任务分配(可用work-pool or processor farm).
1.测试用例
其中底部的数据(real_min, imag_min) to (real_max, imag_max)表示复平面窗口,real_min表示
实部最小值, imag_min表示虚部最小值,real_max表示实部最大值, imag_max表示虚部最大
值.
(-0.84950, 0.21000) to (-0.84860, 0.21090) (0.32000, -0.45000) to (0.50000, 0.05000) (0.26304, 0.00233) to (0.26329, 0.00267)
(-0.63500, 0.68000) to (-0.62500, 0.69000) (-0.46510, -0.56500) to (-0.46470, -0.56460) (-1.50000, -1.00000) to (0.50000, 1.00000)
2.曼德勃罗特集的计算
显示曼德勃罗特集是处理位映射图像的一个例子.首先要对图像进行计算,且计算量很
大.曼德勃罗特集是复数平面中的点集,当对一个函数迭代计算时,这些点将处于准稳定状
态(quasi-stable),即将会增加或减少,但不会超过某一限度.通常该函数为:
zk+1 = zk
2 + c
式中zk+1是复数z = a + bi的第k+1次迭代,c是确定该点在复平面中位置的复数值.z的初始
值为0.迭代将一直进行下去,直到z的幅值(向量长度,这里为22ba+)大于2或者迭代次
数已经达到某种任意的规定的限度.
化简计算.
z2 = a2 - b2 + 2abi
用zreal表示z的实部,zimag表示z的虚部.则
zreal = a2 - b2 + creal
zimag = 2ab + cimag
3.顺序代码
structure complex {
float real;
float image;
};
对一点值的计算并返回迭代次数的
int cal_pixel(complex c)
{ int count, max;
complex z;
float temp, lengthsq;
max = 256;
z.real = 0;
z.imag = 0;
count = 0;
do {
temp = z.real * z.real - z.imag * z.imag + c.real;
z.imag = 2 * z.real * z.imag + c.imag;
z.real = temp;
lengthsqc = z.real * z.real + z.imag * z.imag;
count++;
} while ((lengthsq < 4.0) && (count < max)); //直到z的幅值大于2或者迭代次数达到max
return count;
}
因此,所有给定的曼德勃罗特点必将处在以原点为中心,半径为2的圆中.
计算和显示这些点的代码需要对坐标系统进行一定的缩放来与显示区域的坐标系统相匹配.
假设显示高度为disp_height,宽度为disp_width.而点在显示区域中的位置为(x, y).如果显
示复数平面的这个窗口具有最小值(real_min, imag_min)和最大值(real_max, imag_max),则每
个点需用以下系数加以缩放.
c.real = real_min + x * (real_max - real_min) / disp_width;
c.imag = imag_min + y * (imag_max - imag_min) / disp_height;
设置scale_real = (real_max - real_min) / disp_width;
scale_imag = (imag_max - imag_min) / disp_height;
缩放代码:
for (x = 0; x < disp_width; x++)
for (y = 0; y < disp_height; y++) {
c.real = real_min + ((float) x * scale_real);
c.imag = imag_min + ((float) y * scale_imag);
color = cal_pixel(c);
display(x, y, color);
}
4.并行代码
1)静态任务分配
假定显示区域为640 × 480,在一个进程中要计算10行.即将10 × 640像素变为一组,共48
个进程.为代码如下:
主进程
for ( i = 0, row = 0; i < 48; i++, row = row + 10)
send (&row, pi)
for (i = 0; i < (480 * 640); i++)
recv(&c, &color, PANY);
display(c, color);
}
从进程Slave (process i)
recv(&row, Pmaster);
for (x = 0; x < disp_width; x++)
for (y = row; y < (row + 10); y++) {
c.real = real_min + ((float) x * scale_real);
c.imag = imag_min + ((float) y * scale_imag);
color = cal_pixel(c);
send(&c, &color, Pmaster);
}
改进:成组发送数据.减少通信启动次数.先将结果保存在数组中,然后以一个消息将整个
数组发送给主进程.主进程可用一个通配符以任意顺序接收来自从进程的消息.
2)动态任务分配—工作池/处理器农庄(work pool / processor farm)
动态负载平衡可以用工作池方法实现.在我们的问题中,像素集(更确切应该是坐标集)构成
了任务.任务数是固定的,要计算的像素数在计算前是已知的.各个处理器从工作池中请求
下一个像素子集的坐标.
主进程
count = 0;
row = 0;
for (k = 0; k< procno; k++) {
send(&row, pk, data_tag);
count++;
row++;
}
do {
recv(&slave, &r, color, PANY, result_tag);
count--;
if (row 0)
从进程
recv(y, Pmaster, ANYTAG, source_tag); //接收Pmaster发送的第y行的点
while(source_tag == data_tag) { //判断是否还有消息需要处理
c.imag = imag_min + ((float) y * scale_imag);
for (x = 0; x < disp_width; x++) {
c.real = real_min + ((float) x * scale_real);
color = cal_pixel(c);
}
send(&i, &y, color, Pmaster, result_tag); //将所计算的第y行点的color发给Pmaster
recv(y, &r, Pmaster, source_tag); //接收下一任务
} //如果退出while循环, 则一定是source_tag = termiate_tag, 表明没有任务,程序终止