Mandelbrot 并行实现

最新推荐文章于 2022-01-16 22:45:42 发布

原创最新推荐文章于 2022-01-16 22:45:42 发布 · 5.6k 阅读

13 ·

CC 4.0 BY-SA版权

文章标签：

#pthreads #float #struct #算法 #distance #generator

本文介绍了Mandelbrot集合的概念、定义及与Julia集合的区别，并通过编程实现Mandelbrot集合的并行计算，探讨了使用MPI和pthreads两种方法，展示了并行计算的效果和遇到的问题。

部署运行你感兴趣的模型镜像

最近要交并行计算的作业了，这周终于把作业写了个大概，这期间感觉学了不少东西，总结一下。

Mandelbrot Set 背景

前几天逛维基百科的时候看到了如下的消息：著名数学家、分形之父Benoît B. Mandelbrot（中文名本华·曼德博）美国时间10月15日辞世，享年85岁。

“1979年，在哈佛大学作为访问学者的期间，曼德博开始研究分形集之一——在复平面上一定变换下具有不变性的朱利亚集合。在加斯顿·朱利亚和皮埃尔·法图学术成果的基础上，曼德博利用公式 f_c(z) = z² + c 反复迭代，在计算机上作出了朱利亚集合的图形。在研究朱利亚集合的拓扑结构是怎样依赖于复参数μ的同时，他还提出了后来一他的名字命名的曼德博集合（曼德博集合今天常常由 z² + c 定义，因此曼德博早期以μ为参数所作的图像与后来以 c 为参数所得图像恰如镜面反射般左右对称）。”—From Wikipedia

Wikipedia上Mandelbrot的定义为：

曼德博集合可以用复二次多项式

$f_ c(z) =z^{2}+ c /,$ 来定义

其中c是一个复参数。对于每一个c，从 $z = 0/,$ 开始对f_c(z)进行迭代。

序列 $(0, f_ c(0), f_c(f_ c(0)), f_ c(f_ c(f_ c(0))), /ldots)$ 的值或者延伸到无限大，或者只停留在有限半径的圆盘内。

曼德博集合就是使以上序列不延伸至无限大的所有c点的集合。

公式为： where

从数学上来讲，曼德博集合是一个复数的集合。一个给定的复数c或者属于曼德博集合M，或者不是。

要记住上面介绍的最后一句话：Mandelbrot Set是一个集合，谁的集合？是关于 $f_ c(z) =z^{2}+ c /,$ 当中复数c的集合。刚开始我就老搞不明白Mandelbrot Set和Julia Set的区别，其实他们虽然迭代式大体相同，但是所表示的集合范围却是两码事，下面给出Julia Set的定义：

朱利亚集合可以由下式进行反复迭代得到：

f_c(z) = z² + c

对于固定的复数c，取某一z值（如z = z₀），可以得到序列

$z_0, f_c(z_0), f_c(f_c(z_0)), f_c(f_c(f_c(z_0))), /ldots .$

这一序列可能反散于无穷大或始终处于某一范围之内并收敛于某一值。我们将使其不扩散的z值的集合称为朱利亚集合。

这里需要注意的是：这个迭代式当中，复数c的值是固定的，而Julia的集合的范围是关于z值的范围。

简单来说，对于迭代式：f_c(z) = z² + c ，Mandelbrot Set 迭代过程中的z值是固定的，是使上述迭代式始终在某一范围内而不发散于无穷大的c值的集合；而Julia Set的迭代过程中c是固定的，上述是使上述迭代式始终在某一范围内而不发散于无穷大的z值的集合。

理解了这些定义，就很容易通过编程来实现这样的集合。

最后就是关于绘图着色的问题，大部分都是通过让||z||>R的迭代次数来决定所绘制的颜色，我们定义如果在最大迭代次数N内没有超过R的集合的颜色为黑色，当最大的迭代次数N足够大时，可以认为这些点属于Mandelbrot Set，具体要多大？理论上来说当N为无穷大时，我们绘制的Mandelbrot Set才是真的符合定义，但是我们即使借助于计算机，也不可能处理这种无限的问题。所以说我们定义在最大迭代次数N之内z的模不超过R就认为就是属于Mandelbrot Set其实是不严格的，所幸的是如果z不属于Mandelbrot Set的时候，||z||的增长速度是非常的快，再加上我们一般电脑的屏幕分辨率的限制，当N值足够大的时候所显示的图像与N趋于无穷时显示的图像几乎没有区别。所以N值只要取足够大就可以了。还有一个问题就是R的大小的选择，我查阅了相关的资料，也没仔细看，相关的证明公式太多，懒得看，直接上结论：“If the absolute value of Z (that is, its distance from 0+0i) ever gets bigger than 2, it will never return to a place closer than 2, but it will actually rapidly escape to infinity.”也就是说如果z的模值大于2的时候，可以证明经过上述公式的迭代，z的模值会很快趋于无穷大。所以现在基本都是取R=2。

所以，经过上述的一番理论知识，我们用绘图的方法来展示Mandelbrot Set只是一种演示，并不完全符合Mandelbrot Set的定义，但是对于展示这个集合的相关特性，已经足够了。

刚开始我定义迭代次数n>N时，为黑色，属于Mandelbrot Set，对于不属于Mandelbrot Set的，我还是通过迭代次数n来决定其颜色，当n从N到N/2时，颜色从黑变红，当n从N/2到0时，颜色从红到白，结果出现如下后果，十分血腥。。。

然后我到网上搜索了相关的配色，最后在Rocky Mountain College的K. Stuart Smith 教授那里找到了一份配色的方案，他是将最大迭代次数设为256，应该是和256色有关，但是具体的色彩算法我也不是很清楚，不过他应该是用的一种叫HSB的配色方案，他写的函数就是将HSB的颜色转换为RGB的，具体代码如下：

   1: void setHSBColor (float hue, float saturate, float bright) {

   2:   // when I wrote this, openGL only liked colors specified by RGB values. The

   3:   // mandelbrot routine generated an HSB color, so I wrote this routine to do

   4:   // the conversion. It sure isn't perfect, but it does a respectable job.

   5:   //

   6:   // I expect that part of this work (but the final openGL call) could be

   7:   // pushed back with the mandelbrot color generator for more speedup.

   8:   //

   9:   float red, green, blue;

  10:   float h = (hue * 256) / 60;

  11:   float p = bright * (1 - saturate);

  12:   float q = bright * (1 - saturate * (h - (int)h));

  13:   float t = bright * (1 - saturate * (1 - (h - (int)h)));

14:

  15:   switch ((int)h) {

  16:   case 0:

  17:     red = bright,  green = t,  blue = p;

  18:     break;

  19:   case 1:

  20:     red = q,  green = bright,  blue = p;

  21:     break;

  22:   case 2:

  23:     red = p,  green = bright,  blue = t;

  24:     break;

  25:   case 3:

  26:     red = p,  green = q,  blue = bright;

  27:     break;

  28:   case 4:

  29:     red = t,  green = p,  blue = bright;

  30:     break;

  31:   case 5:

  32:   case 6:

  33:     red = bright,  green = p,  blue = q;

  34:     break;

  35:   }

  36:   glColor3f (red, green, blue);

  37: }

这样，我们连配色方案也有了，核心的算法通过上面的讲述也很清晰了，我的实现如下：

   1: double calcMandelbrotColor(double x,double y)

   2: {

   3:     complex<double> z(0,0);

   4:     complex<double> c(x,y);

   5:     int n=0;

   6:     while(n<256)

   7:     {

   8:         ++n;

   9:         z=z*z+c;

  10:         if(z.real()*z.real()+z.imag()*z.imag()>4) break;

  11:     }

  12:     return (double)n/256.0;

  13: }

这个函数计算出来的是(x,y)点的h值，通过上面的setHSBColor函数中的算法，就可以知道RGB的值,大体代码如下

   1: h=calcMandelbrotColor(x,y);

   2: setHSBColor(h,0.7,1.0-h*h/1.2);//为啥这么用，暂时不知道，应该是和色彩学相关的算法

当然，那个setHSBColor函数中的RGB的值依然是对应于OpenGL中的RGB，每个原色的取值范围为0~1.0,如果用微软的API实现，则要自己转换成0~255,并且自己实现相关的画笔设置功能。

最终的效果如下图：

如果放大后，则有惊人的效果：

既然上面已经提到了Julia Set，那么也顺便实现了一下，其实代码和Mandelbrot Set很相似，区别就是在计算颜色的h值那个地方：

   1: double calcJuliaColor(double x,double y)

   2: {

   3:     complex<double> z(x,y);

   4:     complex<double> c(0.109,0.603);

   5:     int n=0;

   6:     while(n<256)

   7:     {

   8:         ++n;

   9:         z=z*z+c;

  10:         if(z.real()*z.real()+z.imag()*z.imag()>4) break;

  11:     }

  12:     return (double)n/256.0;

  13: }

效果如下，c取<0.109,0.603>

c取<0.314,0.628>

c取<-0.314.-0.628>

c取<-0.814,0.17>

数学之美啊~

Mandelbrot Set的并行实现：

实现并行的方法有很多，这里先说2种--MPI和pthreads，以后用其他方法实现再来补充。

首先说MPI，这个强大的工具已经被很多超级计算机所采用，我没有那么多电脑。。只好拿自己的电脑进行模拟。其实思路很简单，就是假设有n个处理器，拿出一个进行数据的分配工作，称为Master进程(如果真的有n台计算机的话，每个计算机上都会有这个进程，共n个，我是用MPI的mpiexec进行模拟的，也称其为进程吧，其实与真实的n个处理器还是有差别的)，然后剩下的n-1个进程进行数据的操作，称为Slave进程。就用这1个Master进程和n-1个Slave进程进行Mandelbrot Set的计算以及绘图工作，涉及到的通信有：（1）Master将纵坐标切割成n-1个，然后将纵坐标row的信息发送到Slave；(2)每个Slave接受到各自的row的坐标后，然后进行Mandelbrot的颜色的计算，范围为row~row+Height/(n-1),并且将计算的每一个坐标点的颜色值发送给Master;(3)Master接受到Slave发送过来的颜色后，对这个坐标进行绘图工作。

主要代码如下：

   1: //主进程

   2: void calcMandelbrotMaster()

   3: {

   4:     int i,row;

   5:     double pre=0;

   6:     int step=windowY/(nprocs-1);

   7:     int remains=windowY-step*(nprocs-1);

   8:     row=0;

   9:     for(i=1;i<nprocs;++i)

  10:     {

  11:         MPI_Send(&row,1,MPI_INT,i,row_tag,MPI_COMM_WORLD);

  12:         row+=step;

  13:     }

  14:     //暂时认为 (nprocs-1)|windowY

  15:     //if(remains) MPI_Send(&row,1,MPI_INT,nprocs-1,row_tag,MPI_COMM_WORLD);

  16:     for(i=0;i<windowX*windowY;++i)

  17:     {

  18:         MPI_Recv(&currentColor,1,ccolor,MPI_ANY_SOURCE,color_tag,MPI_COMM_WORLD,&status);

  19:         //display(currentColor);

  20:         if(currentColor.c!=pre)

  21:         {

  22:             double b = 1.0 - currentColor.c * currentColor.c / 1.2;

  23:             setHSBColor (currentColor.c, 0.7, b);

  24:         }

  25:         glBegin(GL_POINTS);

  26:         glVertex2i (currentColor.x - windowX/2, windowY/2 - currentColor.y);

  27:         //if(myrank)printf("I am at %d rank/n",myrank);

  28:         glEnd();

  29:         glFlush();

  30:     }

  31: }

   1: //从进程

   2: void calcMandelbrotSlave()

   3: {

   4:     double h,b,pre=0;//pre是前一个颜色，一个小优化

   5:     int row;

   6:     MPI_Recv(&row,1,MPI_INT,0,row_tag,MPI_COMM_WORLD,&status);

   7:     int x,y;

   8:     for(x=0;x<windowX;++x)

   9:         for(y=row;y<min(row+windowY/(nprocs-1),windowY);++y)

  10:         {

  11:             h=calcColor(startX+zoomX*(x),startY+zoomY*(y));//将坐标点传入

  12:             b = 1.0 - h * h / 1.2;

  13:             setHSBColor (h, 0.7, b);

  14:             struct Node cur(x,y,h);

  15:             MPI_Send(&cur,1,ccolor,0,color_tag,MPI_COMM_WORLD);

  16:         }

  17: }

有个问题就是这个程序的绘图使用OpenGL实现，但是绘图完成之后的CPU的使用率依然是100%，不知道怎么解决这个问题，我估计是glutMainLoop()函数的缘故，但是目前不知道怎么解决，暂时先记下来。运行结果如下图：

再来说pthreads，这个其实是用多线程来模拟多核，也是一种并行的方法。

其实用pthreads相对简单，关于其原理就不多说了，大体讲一下方法：

其实用pthreads实现关键在于

   1: pthread_create (pthread_t * tid,

   2:                             const pthread_attr_t * attr,

   3:                             void *(*start) (void *),

   4:                             void *arg);

这个函数，第一个参数是pthreads_t类型的数组，也就是你要建立几个线程，第二个参数是相关属性，一般为NULL，关键是第三个和第四个参数，第三个参数是个函数指针，第四个参数是函数指针指向的函数的参数。最重要的一点：要对绘图的操作进行上锁！主要代码如下：

   1: typedef struct Rect

   2: {

   3:     int xs;

   4:     int xe;

   5:     int ys;

   6:     int ye;

   7: }Rect;

   8: typedef struct Node

   9: {

  10:     HDC *hDC;

  11:     Rect *curRect;

  12: }Node;

  13: Rect myRect[9];//将屏幕分隔成9块

  14: Rect cursorRect;

15:

16:

  17: void *Draw(void *hDCRc)

  18: {

  19:     Rect *curRect=((Node*)hDCRc)->curRect;

  20:     HDC *hDC=((Node*)hDCRc)->hDC;

  21:     int x,y;

  22:     double h,b;

  23:     for(y=curRect->ye;y>=curRect->ys;--y)

  24:     {

  25:         for(x=curRect->xe;x>=curRect->xs;--x)

  26:         {

  27:             pthread_mutex_lock(&drawLock);

  28:             h=calcColor(startX+zoomX*(x),startY+zoomY*(y));//将坐标点传入

  29:             b = 1.0 - h * h / 1.2;

  30:             COLORREF crColor=getHSBColor (h, 0.7, b);

31:

  32:             SetPixel(*hDC,x,y,crColor);

  33:             pthread_mutex_unlock(&drawLock);

  34:         }

  35:     }

  36:     return NULL;

  37: }

  38: void display()

  39: {

  40:     PAINTSTRUCT ps;

  41:     HDC hDC;

  42:     hDC=BeginPaint(hwnd,&ps);

  43:     pthread_t threads[9];

  44:     int i;

  45:     for(i=0;i<9;++i)

  46:     {

  47:         Node *myNode=new Node;

  48:         myNode->curRect=&myRect[i];

  49:         myNode->hDC=&hDC;

  50:         pthread_create(&threads[i],NULL,Draw,myNode);

  51:     }

  52:     for(i=0;i<9;++i)

  53:         pthread_join(threads[i],NULL);

  54:     EndPaint(hwnd,&ps);

  55: }

第一次利用Win32的API进行绘图，遇到了很多问题。比如响应WM_PAINT消息，BeginPaint后一定要EndPaint，要不然会长期占用cpu。暂时记这么多，运行效果如下：

比较有成就感，嘿嘿~

您可能感兴趣的与本文相关的镜像

Anything-LLM

AI应用

AnythingLLM是一个全栈应用程序，可以使用商用或开源的LLM/嵌入器/语义向量数据库模型，帮助用户在本地或云端搭建个性化的聊天机器人系统，且无需复杂设置

Mandelbrot 并行实现

2 条评论