OpenMP 分形生成器

Question

我一直在尝试创建 julia 集的 openMP 变体，但在运行多个线程时无法创建连贯的图像，我一直在尝试解决看起来像竞态但无法解决的问题找出错误。 有问题的输出看起来像所需的输出以及整个图片的“扫描线”。

如果代码不够清楚，我也附上代码。

#include <iostream>
#include <math.h>
#include <fstream>
#include <sstream>
#include <omp.h>
#include <QtWidgets>
#include <QElapsedTimer>
using namespace std;



double newReal(int x, int imageWidth){
    return 1.5*(x - imageWidth / 2)/(0.5 * imageWidth);
}

double newImaginary(int y, int imageHeight){
    return (y - imageHeight / 2) / (0.5 * imageHeight);
}

int julia(double& newReal, double& newImaginary, double& oldReal, double& oldImaginary, double cRe, double cIm,int maxIterations){
    int i;
    for(i = 0; i < maxIterations; i++){
      oldReal = newReal;
      oldImaginary = newImaginary;
      newReal = oldReal * oldReal - oldImaginary * oldImaginary + cRe;
      newImaginary = 2 * oldReal * oldImaginary + cIm;
      if((newReal * newReal + newImaginary * newImaginary) > 4) break;
    }
    return i;
}
int main(int argc, char *argv[])
{
    int fnum=atoi(argv[1]);
    int numThr=atoi(argv[2]);
//    int imageHeight=atoi(argv[3]);
//    int imageWidth=atoi(arg[4]);
//    int maxIterations=atoi(argv[5]);
//    double cRe=atof(argv[3]);
//    double cIm=atof(argv[4]);
    //double cRe, cIm;
    int imageWidth=10000, imageHeight=10000, maxIterations=3000;
    double newRe, newIm, oldRe, oldIm,cRe,cIm;
    cRe = -0.7;
    cIm = 0.27015;
    string fname;
    QElapsedTimer time;
    QImage img(imageHeight, imageWidth, QImage::Format_RGB888);//Qimagetesting
    img.fill(QColor(Qt::black).rgb());//Qimagetesting
    time.start();
    int i,x,y;
    int r, gr, b;
#pragma omp parallel for shared(imageHeight,imageWidth,newRe,newIm) private(x,y,i) num_threads(3)
        for(y = 0; y < imageHeight; y++)
        {
            for(x = 0; x < imageWidth; x++)
            {
                newRe = newReal(x,imageWidth);
                newIm = newImaginary(y,imageHeight);
                i= julia(newRe, newIm, oldRe, oldIm, cRe, cIm, maxIterations);

                r = (3*i % 256);
                gr = (2*(int)sqrt(i) % 256);
                b = (i % 256);
                img.setPixel(x, y, qRgb(r, gr, b));
             }
        }

    //stringstream s;
    //s << fnum;
    //fname= "julia" + s.str();
    //fname+=".png";
    //img.save(fname.c_str(),"PNG", 100);
    img.save("julia.png","PNG", 100);
    cout<< "Finished"<<endl;
    cout<<time.elapsed()/1000.00<<" seconds"<<endl;
}

Answer 1

正如 cmets 所指出的，您有两个主要问题：

• newRe 和 newIm 是共享的，但不应该是
• 未指定r、gr 和b 的访问权限（我认为默认共享）
• 有并发调用QImage::setPixel

要纠正此问题，请不要犹豫，将 omp for 循环嵌套在 omp parallel 块中。

在 for 循环之前声明私有变量：

为防止同时调用QImage::setPixel、since this function is not thread safe，可以将其置于临界区，使用#pragma omp critical。

---

int main(int argc, char *argv[])
{
    int imageWidth=1000, imageHeight=1000, maxIterations=3000;

    double cRe = -0.7;
    double cIm = 0.27015;

    QElapsedTimer time;
    QImage img(imageHeight, imageWidth, QImage::Format_RGB888);//Qimagetesting
    img.fill(Qt::black);
    time.start();

    #pragma omp parallel
    {
        /* all folowing values will be private */
        int i,x,y;
        int r, gr, b;
        double newRe, newIm, oldRe, oldIm;

        #pragma omp for
        for(y = 0; y < imageHeight; y++)
        {            
            for(x = 0; x < imageWidth; x++)
            {
                newRe = newReal(x,imageWidth);
                newIm = newImaginary(y,imageHeight);
                i= julia(newRe, newIm, oldRe, oldIm, cRe, cIm, maxIterations);

                r = (3*i % 256);
                gr = (2*(int)sqrtf(i) % 256);
                b = (i % 256);

                #pragma omp critical
                    img.setPixel(x, y, qRgb(r, gr, b));
             }
        }
    }

    img.save("julia.png","PNG", 100);

    cout<<time.elapsed()/1000.00<<" seconds"<<endl;

    return 0;
}

---

更进一步，您可以节省一些 CPU 时间，将 ::setPixel 替换为 ::scanLine：

    #pragma omp for
    for(y = 0; y < imageHeight; y++)
    {                        
        uchar *line = img.scanLine(y);
        for(x = 0; x < imageWidth; x++)
        {
            newRe = newReal(x,imageWidth);
            newIm = newImaginary(y,imageHeight);
            i= julia(newRe, newIm, oldRe, oldIm, cRe, cIm, maxIterations);

            r = (3*i % 256);
            gr = (2*(int)sqrtf(i) % 256);
            b = (i % 256);

            *line++ = r;
            *line++ = gr;
            *line++ = b;
         }
    }

---

编辑：

由于 julia 集似乎在(0,0) 点周围有一个中心对称性，因此您只能执行一半的微积分：

int half_heigt = imageHeight / 2;
#pragma omp for

// compute only for first half of image
for(y = 0; y < half_heigt; y++)
{            
    for(x = 0; x < imageWidth; x++)
    {
        newRe = newReal(x,imageWidth);
        newIm = newImaginary(y,imageHeight);
        i= julia(newRe, newIm, oldRe, oldIm, cRe, cIm, maxIterations);

        r = (3*i % 256);
        gr = (2*(int)sqrtf(i) % 256);
        b = (i % 256);

        #pragma omp critical
        {
            // set the point  
            img.setPixel(x, y, qRgb(r, gr, b));

            // set the symetric point
            img.setPixel(imageWidth-1-x, imageHeight-1-y, qRgb(r, gr, b));
        }
     }
}