我が家は週末になるとセロハンテープがなくなります。主に次男がやたらと工作した車や剣などにベタベタ張るのが原因です。長男は無駄遣いすることはほとんどないようですが、大好きなスポーツカーをたくさんつくるので、消費量は多いです。@kjunichiです。

背景

去年Juliaで並列処理が意外と簡単にできることがわかり、aobenchのjulia版を作り、

aobenchをjuliaでやってみた

という記事を書いて、並列処理もその続きに書こうとはしていたものの、すっかり忘れていたのと、昨年末自作PCをXeonの8コアにしたので、マルチコアの恩恵がより受けられる環境が身近になったのもあり記事にした。

簡単な並列化

単純にJuliaに用意されている並列処理のキーワードを追加して、並列処理させる。

addprocsでプロセス数を設定する。

addprocs(4)

複数のプロセスから参照する変数に@everywhereキーワードをつける。

@everywhere hoge

非同期処理は@async、@syncキーワードで囲むことでJavaなどのjoin的な動きになる。

@sync begin
 @async remotecall_wait(p, render, fimg,img, WIDTH, HEIGHT, NSUBSAMPLES)
end

成果物

# original code:
# https://github.com/syoyo/aobench/blob/master/ao.c
addprocs(7)

@everywhere WIDTH        =512
@everywhere HEIGHT       =512
@everywhere NSUBSAMPLES  =2
@everywhere NAO_SAMPLES =8

@everywhere type Vec
    x::Float64
    y::Float64
    z::Float64
end

@everywhere type Isect
    t::Float64
    p::Vec
    n::Vec
    hit::Int64
end

@everywhere type Sphere
    center::Vec
    radius::Float64
end

@everywhere type Plane
    p::Vec
    n::Vec
end

@everywhere type Ray
    org::Vec
    dir::Vec
end

@everywhere function vdot(v0, v1)
    v0.x * v1.x + v0.y * v1.y + v0.z * v1.z
end

@everywhere function vcross(v0, v1)
    Vec(v0.y * v1.z - v0.z * v1.y,
               v0.z * v1.x - v0.x * v1.z,
                 v0.x * v1.y - v0.y * v1.x)
end

@everywhere function vnormalize(c)

    length = sqrt(vdot(c, c))

    if (abs(length) > 1.0e-17)
        c.x /= length
        c.y /= length
        c.z /= length
    end
end

@everywhere function ray_sphere_intersect(isect, ray, sphere)

    rs = Vec(ray.org.x - sphere.center.x,
                     ray.org.y - sphere.center.y,
                       ray.org.z - sphere.center.z)

    B = vdot(rs, ray.dir)
    C = vdot(rs, rs) - sphere.radius * sphere.radius
    D = B * B - C

    if (D > 0.0)
        t = -B - sqrt(D)

        if ((t > 0.0) && (t < isect.t))
            isect.t = t
            isect.hit = 1

            isect.p.x = ray.org.x + ray.dir.x * t
            isect.p.y = ray.org.y + ray.dir.y * t
            isect.p.z = ray.org.z + ray.dir.z * t

            isect.n.x = isect.p.x - sphere.center.x
            isect.n.y = isect.p.y - sphere.center.y
            isect.n.z = isect.p.z - sphere.center.z

            vnormalize(isect.n)
        end
    end
end

@everywhere function ray_plane_intersect(isect, ray, plane)

    d = -vdot(plane.p, plane.n)
    v = vdot(ray.dir, plane.n)

    if (abs(v) < 1.0e-17)
     return
    end

    t = -(vdot(ray.org, plane.n) + d) / v;

    if ((t > 0.0) && (t < isect.t))
        isect.t = t
        isect.hit = 1

        isect.p.x = ray.org.x + ray.dir.x * t
        isect.p.y = ray.org.y + ray.dir.y * t
        isect.p.z = ray.org.z + ray.dir.z * t

        isect.n = plane.n
    end
end

@everywhere function orthoBasis(basis, n)

    basis[3] = n
    basis[2].x = 0.0; basis[2].y = 0.0; basis[2].z = 0.0

    if (n.x < 0.6) && (n.x > -0.6)
        basis[2].x = 1.0
    elseif ((n.y < 0.6) && (n.y > -0.6))
        basis[2].y = 1.0
    elseif ((n.z < 0.6) && (n.z > -0.6))
        basis[2].z = 1.0
    else
        basis[2].x = 1.0
    end

    basis[1] = vcross(basis[2], basis[3])
    vnormalize(basis[1])

    basis[2] = vcross(basis[3], basis[1])
    vnormalize(basis[2])
end

@everywhere function drand48()
  rand()
end

@everywhere function ambient_occlusion(isect)

    ntheta = NAO_SAMPLES
    nphi   = NAO_SAMPLES
    eps = 0.0001;

    p =Vec(
             isect.p.x + eps * isect.n.x,
               isect.p.y + eps * isect.n.y,
                     isect.p.z + eps * isect.n.z)

    basis = [Vec(0.0,0.0,0.0),Vec(0.0,0.0,0.0),Vec(0.0,0.0,0.0)]
    orthoBasis(basis, isect.n)

    occlusion = 0.0

    for j = 0:ntheta-1
        for i = 0:nphi-1
            theta = sqrt(drand48())
            phi   = 2.0 * pi * drand48()

            x = cos(phi) * theta
            y = sin(phi) * theta
            z = sqrt(1.0 - theta * theta)

            # local . global
             rx = x * basis[1].x + y * basis[2].x + z * basis[3].x
             ry = x * basis[1].y + y * basis[2].y + z * basis[3].y
             rz = x * basis[1].z + y * basis[2].z + z * basis[3].z

            ray = Ray(p,Vec(rx,ry,rz))

            occIsect = Isect(1.0e+17,Vec(0.0,0.0,0.0),Vec(0.0,0.0,0.0),0)

            ray_sphere_intersect(occIsect, ray, spheres[1])
            ray_sphere_intersect(occIsect, ray, spheres[2])
            ray_sphere_intersect(occIsect, ray, spheres[3])
            ray_plane_intersect(occIsect, ray, plane)

            if occIsect.hit != 0;
             occlusion += 1.0
            end

        end
    end

    occlusion = (ntheta * nphi - occlusion) / Float64(ntheta * nphi)
    Vec(occlusion, occlusion, occlusion)
end

@everywhere function clamp(f::Float64)
  i = Int64(floor(f * 255.5))

  if (i < 0) i = 0;end
  if (i > 255) i = 255;end
  Cuchar(i);
end

@everywhere function render(fimg::SharedArray,img, w, h, nsubsamples)

  # find startx,starty
  idx = indexpids(fimg)
  @show idx
  if idx == 0
    startx = 0
    starty = 0
  else
    nchunks = length(procs(fimg))
    unitsize = round(Int, WIDTH*HEIGHT*3/nchunks)
    startpos = unitsize * (idx-1)
    @show nchunks, unitsize, startpos
    starty = floor(Int,startpos/(WIDTH*3))
    startx = startpos % WIDTH
  end
  @show startx, starty
  iidx = 0
  for y = starty:h-1
      for x = startx:w-1

          for v = 0:nsubsamples-1
              for u = 0:nsubsamples-1
                  px = (x + (u / Float64(nsubsamples)) - (w / 2.0)) / (w / 2.0);
                  py = -(y + (v / Float64(nsubsamples)) - (h / 2.0)) / (h / 2.0);

                  ray = Ray(Vec(0.0,0.0,0.0),Vec(px,py,-1.0))

                  vnormalize(ray.dir)

                  isect = Isect(1.0e+17,Vec(0.0,0.0,0.0),Vec(0.0,0.0,0.0),0)

                  ray_sphere_intersect(isect, ray, spheres[1]);
                  ray_sphere_intersect(isect, ray, spheres[2]);
                  ray_sphere_intersect(isect, ray, spheres[3]);
                  ray_plane_intersect(isect, ray, plane);

                  if isect.hit != 0
                      col = ambient_occlusion(isect);

                      fimg[3 * (y * w + x) + 1] += col.x;
                      fimg[3 * (y * w + x) + 2] += col.y;
                      fimg[3 * (y * w + x) + 3] += col.z;
                  end
              end
          end

          fimg[3 * (y * w + x) + 1] /= Float64(nsubsamples * nsubsamples);
          fimg[3 * (y * w + x) + 2] /= Float64(nsubsamples * nsubsamples);
          fimg[3 * (y * w + x) + 3] /= Float64(nsubsamples * nsubsamples);

          img[3 * (y * w + x) + 1] = clamp(fimg[3 *(y * w + x) + 1]);
          img[3 * (y * w + x) + 2] = clamp(fimg[3 *(y * w + x) + 2]);
          img[3 * (y * w + x) + 3] = clamp(fimg[3 *(y * w + x) + 3]);
          iidx+=3
          if(iidx>unitsize)
            return
          end
      end
  end
end

@everywhere function saveppm(fname,  w,  h, img)
    open(fname, "w") do fp
      write(fp, "P6\n");
      write(fp, "$w $h\n");
      write(fp, "255\n");
      write(fp, img)
    end
end

fimg = SharedArray(Float64,WIDTH * HEIGHT * 3);
zeros(fimg)

img = SharedArray(Cuchar,WIDTH * HEIGHT * 3)

# init_scene
@everywhere spheres =[
                Sphere(Vec(-2.0,0.0,-3.5),0.5),
              Sphere(Vec(-0.5,0.0,-3.0),0.5),
                Sphere(Vec(1.0,0.0,-2.2),0.5)]
@everywhere plane = Plane(Vec(0.0,-0.5,0.0),Vec(0.0,1.0,0.0))
@sync begin
  for p in procs(fimg)
    @async remotecall_wait(render,p, fimg,img, WIDTH, HEIGHT, NSUBSAMPLES)
  end
end
#render(fimg,img, WIDTH, HEIGHT, NSUBSAMPLES);

saveppm("ao.ppm", WIDTH, HEIGHT, img);

ほんとに速くなったのか

以下、ao.jlがシングルコア版、aop.jlがマルチプロセス版

C:\Users\kjw_j\Documents\work\julia
λ  Measure-Command {julia .\ao.jl}


Days              : 0
Hours             : 0
Minutes           : 0
Seconds           : 26
Milliseconds      : 479
Ticks             : 264795917
TotalDays         : 0.000306476755787037
TotalHours        : 0.00735544213888889
TotalMinutes      : 0.441326528333333
TotalSeconds      : 26.4795917
TotalMilliseconds : 26479.5917



C:\Users\kjw_j\Documents\work\julia
λ  Measure-Command {julia .\aop.jl}


Days              : 0
Hours             : 0
Minutes           : 0
Seconds           : 28
Milliseconds      : 862
Ticks             : 288628742
TotalDays         : 0.000334061043981481
TotalHours        : 0.00801746505555555
TotalMinutes      : 0.481047903333333
TotalSeconds      : 28.8628742
TotalMilliseconds : 28862.8742



C:\Users\kjw_j\Documents\work\julia
λ

マルチプロセス版の方がおそい！