2013-04-21 2 views
1

У меня есть следующий заданный код Мандельброта в OpenMP. Мой код C работает очень хорошо, и изображение, которое оно производит, идеально. Но используя OpenMP, он компилируется и работает правильно, но, к сожалению, я не могу открыть выходной файл .ppm, просто Gimp не может его прочитать.Нужна корректировка моего установленного кода OpenMP Mandelbrot

// mandopenmp.c 
// to compile: gcc -fopenmp mandopenmp.c -o mandopenmp -lm 
// usage: ./mandopenmp <no_of_iterations> > output.ppm 

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <time.h> 
#include <omp.h> 

typedef struct { 
    int r, g, b; 
} rgb; 


void color(rgb **m, int x, int y, int red, int green, int blue) 
{ 
    m[x][y].r = red; 
    m[x][y].g = green; 
    m[x][y].b = blue; 
} 

void mandelbrot(int niterations, rgb **m) 
{ 
    int w = 600, h = 400, x, y, i; 
    // each iteration, it calculates: newz = oldz*oldz + p, 
    // where p is the current pixel, and oldz stars at the origin 
    double pr, pi;     // real and imaginary part of the pixel p 
    double newRe, newIm, oldRe, oldIm; // real and imaginary parts of new and old z 
    double zoom = 1, moveX = -0.5, moveY = 0; // you can change these to zoom and change position 

    printf("P6\n# AUTHOR: Erkan Tairi\n"); 
    printf("%d %d\n255\n",w,h); 

    //loop through every pixel 
    #pragma omp parallel for private(x,i,pr,pi,newRe,newIm,oldRe,oldIm) schedule(dynamic, 1) 
    for(y = 0; y < h; y++) { 
     for(x = 0; x < w; x++) { 
      // calculate the initial real and imaginary part of z, 
      // based on the pixel location and zoom and position values 
      pr = 1.5 * (x - w/2)/(0.5 * zoom * w) + moveX; 
       pi = (y - h/2)/(0.5 * zoom * h) + moveY; 
       newRe = newIm = oldRe = oldIm = 0; //these should start at 0,0 
       // start the iteration process 
       for(i = 0; i < niterations; i++) { 
         // remember value of previous iteration 
         oldRe = newRe; 
         oldIm = newIm; 
         // the actual iteration, the real and imaginary part are calculated 
         newRe = oldRe * oldRe - oldIm * oldIm + pr; 
         newIm = 2 * oldRe * oldIm + pi; 
         // if the point is outside the circle with radius 2: stop 
         if((newRe * newRe + newIm * newIm) > 4) break; 
       } 
       if(i == niterations) 
       color(m, x, y, 0, 0, 0); // black 
      else 
      { 
       // normalized iteration count method for proper coloring 
       double z = sqrt(newRe * newRe + newIm * newIm); 
       int brightness = 256. * log2(1.75 + i - log2(log2(z)))/log2((double)niterations); 
       color(m, x, y, brightness, brightness, 255); 
      } 
      } 
    } 
} 

int main(int argc, char *argv[]) 
{ 
    int niterations, i, j; 

    if(argc != 2) 
    { 
     printf("Usage: %s <no_of_iterations> > output.ppm\n", argv[0]); 
     exit(1); 
    } 

    niterations = atoi(argv[1]); 

    rgb **m; 
    m = malloc(600 * sizeof(rgb *)); 
    for(i = 0; i < 600; i++) 
     m[i] = malloc(400 * sizeof(rgb)); 

    double begin = omp_get_wtime(); 
    mandelbrot(niterations, m); 

    for(i = 0; i < 600; i++) { 
     for(j = 0; j < 400; j++) { 
      fputc((char)m[i][j].r, stdout); 
      fputc((char)m[i][j].g, stdout); 
      fputc((char)m[i][j].b, stdout); 
     } 
    } 

    double end = omp_get_wtime(); 

    double time_spent = end - begin; 
    fprintf(stderr, "Elapsed time: %.2lf seconds.\n", time_spent); 

    for(i = 0; i < 600; i++) 
     free(m[i]); 
    free(m); 

    return 0; 
} 
+0

Вы также можете быть заинтересованы в коде здесь https://stackoverflow.com/questions/48069990/multithreaded-simd-vectorized-mandelbrot-in-r-using-rcpp-openmp –

ответ

4

Я не знаю внутренних элементов набора Mandrelbot, но я сделаю снимок, основанный на рабочем процессе вашей программы.

Возможно, это потому, что вы пишете цвет в выходной файл в параллельном разделе. Это означает, что ваши пиксели записываются по завершении вычислительного процесса, но это не означает, что вычислительный процесс пикселя X закончится до обработки пикселя X+1.

Таким образом, при записи в файл, вы в конечном итоге пишете сначала (например) пиксель X+1, а затем пиксель X, смешивая цвета.

Попробуйте записать результат вывода в матрицу. Вам нужно будет изменить свою функцию color, добавив два параметра: i и j с координатами пикселя для записи.

После того, как вся обработка завершается, и каждый пиксель вычисляется, тогда он должен записать пиксели матрицы в выходной файл.

Код:

typedef struct { 
    int r, g, b; 
} rgb; 

void color(rgb **m, int x, int y, int red, int green, int blue) { 
    m[x][y].r = red; 
    m[x][y].g = green; 
    m[x][y].b = blue; 
} 

void mandelbrot(rgb **m, int niterations) { // note the new argument, m. 
    // and your code goes on and on... until: 
      if (i == niterations) 
       color(m, x, y, 0, 0, 0); 
      else { 
       // normalized iteration count method for proper coloring 
       double z = sqrt(newRe * newRe + newIm * newIm); 
       int brightness = 256. * log2(1.75 + i - log2(log2(z)))/log2((double)niterations); 
       color(m, x, y, brightness, brightness, 255); 
      } 
     } 
    } 
} 

int main(int argc, char *argv[]) { 
    // everything ok until... 

    double begin = omp_get_wtime(); 

    rgb **m; 
    m = malloc(sizeof(rgb*) * 600); 
    for (i = 0; i < 600; i++) { 
     m[i] = malloc(400 * sizeof(rgb)); 

    // finally call mandelbrot! 
    mandelbrot(m, niterations); 
    double end = omp_get_wtime(); 

    // now that you have computed your set, you just walk the array writing the output to the file. 

    for (i = 0; i < 600; i++) { 
     free(m[i]); 
    } 
    free(m); 

    double time_spent = end - begin; 
    fprintf(stderr, "Elapsed time: %.2lf seconds.\n", time_spent); 

    return 0; 
} 
+0

Да, я думал это, в то время как параллельно любой поток может завершить свое выполнение перед другим потоком. Но я действительно не знаю, как это исправить. Возможно ли вам, написать короткий пример кода о том, что вы имеете в виду? – 2013-04-21 22:31:20

+0

Отправленный код. Не очень хорошо написано, но идея дана (извините, я тороплюсь ... но мне все равно нужна помощь :)). Пожалуйста, прокомментируйте здесь, если у вас остались вопросы. –

+0

Какой будет размер моей матрицы, это будет мой размер изображения, то есть 600x400? Я редактировал оригинальный код и делал некоторые глупые изменения. Хотя теперь я получаю много ошибок. Я предполагаю, что они в основном касаются выделения и освобождения моей матрицы. – 2013-04-21 23:06:37

0

Если вы пишете в файл последовательно (что вы делали в исходном коде, прежде чем редактировать его), то вы можете использовать ordered прагму, прежде чем писать в файл , Это сделает изображение правильным с использованием исходного кода. См. Следующую ссылку http://bisqwit.iki.fi/story/howto/openmp/#ExampleCalculatingTheMandelbrotFractalInParallel

Однако это не оптимальное решение. Оптимальные решения - сначала записать в буфер памяти, а после того, как код mandelbrot завершит заполнение буфера, напишите в файл (как и в новом коде).

У меня есть несколько предложений по ускорению вашего кода. Предохраните цикл x и y (как показано в ссылке) и используйте динамическое расписание (также показано в этой ссылке), поскольку каждый пиксель занимает разное время. Наконец, используйте SSE/AVX для работы на двух (SSE) или четырех (AVX) пикселах одновременно. В общей сложности вы должны получить более 20 раз с OpenMP и SSE/AVX.

#pragma omp ordered { 
    if(i == niterations) 
     color(m, x, y, 0, 0, 0); // black - use original color function which writes to file 
    else 
    { 
     // normalized iteration count method for proper coloring 
     double z = sqrt(newRe * newRe + newIm * newIm); 
     int brightness = 256. * log2(1.75 + i - log2(log2(z)))/log2((double)niterations); 
     color(m, x, y, brightness, brightness, 255); //use original color function which writes to file 
    } 
} 
+0

Я видел ссылку, спасибо за это. Можете ли вы помочь мне добиться большего ускорения? – 2013-04-23 16:30:08

+0

Я проверил бы следующую ссылку http://www.bealto.com/mp-mandelbrot.html. Однако он не использует AVX, так что вы можете немного ускорить его в CPU. У меня есть код, который делает это, но я бы сначала обработал ваш код масштабирования. Я планирую разместить сайт с моими результатами в какой-то момент ... – 2013-04-23 19:39:51

1

Ваша реализация испорчена. Вы указали много переменных, которые должны быть private вместо shared. Сюда входят pr, pi, newRe, newIm. Также oldRe и oldIm делятся по умолчанию, поскольку они объявлены в области, которая является внешней для параллельной области. Все они должны быть частными вместо:

#pragma omp parallel for private(x,i,pr,pi,newRe,newIm,oldRe,oldIm) 

Также планирование по умолчанию для parallel for петель часто (но не всегда) static. Это не является оптимальным для таких вещей, как фракталы, поскольку для вычисления каждой строки или столбца на изображении требуется различное время. Поэтому вы должны применить предложение schedule(dynamic,1) и играть с размером куска (1 в этом случае), пока не получите лучшее ускорение.

#pragma omp parallel for private(x,i,pr,pi,newRe,newIm,oldRe,oldIm) \ 
      schedule(dynamic,1) 
+0

Спасибо за эту информацию, это помогло мне добиться ускорения. – 2013-04-22 19:47:40