使用OpenMP編譯smallpt導致運行時無限循環

Question

我目前正在查看Keavin Beason的smallpt代碼。我使用g++ -O3 -fopenmp smallpt.cpp編譯了代碼，並且我遇到了似乎是無限循環或死鎖的問題。

僅使用g++ -O3 smallpt.cpp編譯代碼就會在其頁面上看到圖像，但我無法獲得OpenMP並行化的功能。

僅供參考，我使用Cygwin和GCC 4.5.0在Windows 7 64位機器上進行編譯。作者自己已經聲明他運行的是相同的確切代碼，並且沒有遇到任何問題，但是我無法讓程序在完成跟蹤圖像時退出。

這可能是我的特定編譯器和環境的問題，還是我在這裏做錯了什麼？以下是使用OpenMP進行並行化的特定代碼片段。我只修改了一些小的格式以使其更具可讀性。

---

int main(int argc, char *argv[]) 
{ 
    int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1; 

    Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir 
    Vec cx=Vec(w*.5135/h); 
    Vec cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h]; 

    #pragma omp parallel for schedule(dynamic, 1) private(r)  // OpenMP 
    for (int y=0; y<h; y++)      // Loop over image rows 
    { 
    fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1)); 
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++) // Loop cols 
    { 
     for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)  // 2x2 subpixel rows 
     { 
     for (int sx=0; sx<2; sx++, r=Vec())  // 2x2 subpixel cols 
     { 
      for (int s=0; s<samps; s++) 
      { 
      double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1); 
      double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2); 
      Vec d = cx*(((sx+.5 + dx)/2 + x)/w - .5) + 
        cy*(((sy+.5 + dy)/2 + y)/h - .5) + cam.d; 
      r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps); 
      } // Camera rays are pushed ^^^^^ forward to start in interior 
      c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25; 
     } 
     } 
    } 
    } 

    /* PROBLEM HERE! 
     The code never seems to reach here 
     PROBLEM HERE! 
    */ 
    FILE *f = fopen("image.ppm", "w");   // Write image to PPM file. 
    fprintf(f, "P3\n%d %d\n%d\n", w, h, 255); 
    for (int i=0; i<w*h; i++) 
    fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z)); 
} 

下面是輸出的程序產生，當它運行到完成：

$ time ./a 
Rendering (4 spp) 100.00%spp) spp) 00..0026%% 

---

以下是最基本的代碼，可以重現上述行爲

#include <cstdio> 
#include <cstdlib> 
#include <cmath> 

struct Vector 
{ 
    double x, y, z; 
    Vector() : x(0), y(0), z(0) {} 
}; 

int toInt(double x) 
{ 
    return (int)(255 * x); 
} 

double clamp(double x) 
{ 
    if (x < 0) return 0; 
    if (x > 1) return 1; 
    return x; 
} 

int main(int argc, char *argv[]) 
{ 
    int w = 1024; 
    int h = 768; 
    int samples = 1; 

    Vector r, *c = new Vector[w * h]; 

    #pragma omp parallel for schedule(dynamic, 1) private(r) 
    for (int y = 0; y < h; y++) 
    { 
    fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samples * 4, 100. * y/(h - 1)); 
    for (unsigned short x = 0, Xi[3]= {0, 0, y*y*y}; x < w; x++) 
    { 
     for (int sy = 0, i = (h - y - 1) * w + x; sy < 2; sy++) 
     { 
     for (int sx = 0; sx < 2; sx++, r = Vector()) 
     { 
      for (int s = 0; s < samples; s++) 
      { 
      double r1 = 2 * erand48(Xi), dx = r1 < 1 ? sqrt(r1) - 1 : 1 - sqrt(2 - r1); 
      double r2 = 2 * erand48(Xi), dy = r2 < 1 ? sqrt(r2) - 1 : 1 - sqrt(2 - r2); 
      r.x += r1; 
      r.y += r2; 
      } 

      c[i].x += clamp(r.x)/4; 
      c[i].y += clamp(r.y)/4; 
     } 
     } 
    } 
    } 

    FILE *f = fopen("image.ppm", "w");   // Write image to PPM file. 
    fprintf(f, "P3\n%d %d\n%d\n", w, h, 255); 
    for (int i=0; i<w*h; i++) 
    fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z)); 
} 

這是從下面的示例程序獲得的輸出：

$ g++ test.cpp 
$ ./a 
Rendering (4 spp) 100.00% 

$ g++ test.cpp -fopenmp 
$ ./a 
Rendering (4 spp) 100.00%spp) spp) 00..0052%% 

Answer 1

fprintf不受臨界區或#pragma omp single/master保護。如果在Windows上這個東西弄亂了控制檯，我不會感到驚訝。