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

Mandelbrot Set 背景

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

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

Wikipedia上Mandelbrot的定义为：

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

来定义

其中c是一个复参数。对于每一个c，从开始对fc(z)进行迭代。

序列 的值或者延伸到无限大，或者只停留在有限半径的圆盘内。

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

公式为： where

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

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

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

fc(z) = z2 + c

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

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

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

简单来说，对于迭代式：fc(z) = z2 + 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。暂时记这么多，运行效果如下：

比较有成就感，嘿嘿~

Mandelbrot 并行实现相关推荐

  1. 并行程序设计报告(MPI并行计算π，实现mandelbrot集）

    代码的github地址 一.熟悉MPI并行程序设计环境 1.硬件 电脑:HP暗夜精灵 内存:4G 处理器:ntel® Core™ i5-6300HQ CPU @ 2.30GHz × 4 显卡:NVID ...

  2. OpenMP并行化实例----Mandelbrot集合并行化计算

    在理想情况下,编译器使用自动并行化能够管理一切事务,使用OpenMP指令的一个优点是将并行性和算法分离,阅读代码时候无需考虑并行化是如何实现的.当然for循环是可以并行化处理的天然材料,满足一些约束的 ...

  3. 《深入理解并行编程》中文版

    原文的下载地址:http://kernel.org/pub/linux/kernel/people/paulmck/perfbook/perfbook.html 中文版下载地址:深入理解并行编程V1. ...

  4. mpi并行 写同一文件_并行计算调度策略的笔记(001)

    Load Balancing & Termination 因为并行程序的运行速度主要取决于最慢的那个进程.所以保证每个进程的运行时间差不多是非常重要的,这也是负载平衡. 并行计算中的任务调度 ...

  5. Racket编程指南——20 并行

    20 并行 Racket提供两种形式的并行(parallelism):前程(futures)和现场(places).在提供多个处理器的平台上,并行可以提高程序的运行时性能. 关于Racket里顺序性能 ...

  6. 并发 vs 并行 （Concurrency Is Not Parallelism）

    前言 不知你是否曾经下列这些疑问? 并发与并行性有何关系? 什么是同步和异步执行? 如何区分并发与并行? 线程如何与所有这些概念一起使用? 并发 并发性意味着应用程序同时(并发地)在多个任务上取得进展 ...

  7. 分布式训练使用手册-paddle 数据并行

    分布式训练使用手册¶ 分布式训练基本思想¶ 分布式深度学习训练通常分为两种并行化方法:数据并行,模型并行,参考下图: 在模型并行方式下,模型的层和参数将被分布在多个节点上,模型在一个mini-batc ...

  8. 深度学习的分布式训练--数据并行和模型并行

    <div class="htmledit_views"> 在深度学习这一领域经常涉及到模型的分布式训练(包括一机多GPU的情况).我自己在刚刚接触到一机多卡,或者分布式 ...

  9. PyTorch Data Parrallel数据并行

    PyTorch Data Parrallel数据并行 • 可选择:数据并行处理 • 本文将学习如何用 DataParallel 来使用多 GPU. 通过 PyTorch 使用多个 GPU 非常简单.可 ...

最新文章

  1. R语言构建随机森林模型错误解决：Error in y - ymean : non-numeric argument to binary operator
  2. 华为今年不会发布鸿蒙系统的手机,华为：今年不会推出鸿蒙系统手机 将坚守安卓生态...
  3. Python 安装cx_Oracle模块折腾笔记
  4. [python] 3 、基于串口通信的嵌入式设备上位机自动测试程序框架（简陋框架）...
  5. 跨平台的 .NET 运行环境 Mono 3.2 新特性
  6. python学习笔记（二）— 集合
  7. class没有发布到tomcat_Java 类在 Tomcat 中是如何加载的？
  8. boost::mp11::mp_same相关用法的测试程序
  9. iap java md5_苹果应用内支付(iOS IAP)的流程与常用攻击方式
  10. c语言画爱心附带解释,用C语言画一个“爱心”
  11. datagridview 纵向 横向 合并单元格_Excel横向（行）筛选技巧分享，别人3分钟，你只要10秒...
  12. QueryRunner类 的應用，以及ResultSetHandler 接口的实现类
  13. 双层PDFmaker
  14. 数据结构-六度空间（模拟六度分隔理论）
  15. Modulo Sum
  16. 腾讯云-产品开通和密钥查看
  17. systemctl 是管制服务的主要工具
  18. asp.net gridview itemtemplate中控件事件获取行参数
  19. 青花瓷音乐的单片机c语言程序,c语言曲谱_单片机c语言音乐简谱代码
  20. MySQL--事务回滚机制与原理

热门文章

  1. 【VB与数据库】机房收费系统开发阶段之登陆
  2. opa example
  3. Java笔记整理——包装类
  4. python怎么运行代码-python如何运行代码
  5. 技术科普丨马达测试知多少？
  6. 借贷期末余额 oracle,期初余额和期末余额的借贷方是一样的吗
  7. KEIL编程中常见的警告
  8. crc-16 ccitt标准在哪里可以看到详细的
  9. java实施工程师自我介绍
  10. Hadoop MapReduce 统计汽车销售信息