IPython + CTypes + Intel Visual Fortran, Mandelbrot set

Fortran の 配列を Python 側と DLL で共有して引き渡しています。 Mandelbrot 集合の計算は一瞬で終わりますが、Python で配列をタプルの配列に変換するループが遅すぎてこっちの方で地獄のように時間がかかります。I/O をかませたとしても Fortran 側で BMP で出力したものを Python 側で読み込んだ方がよほど早いです。この辺は配列のデータ構造の変換の問題なので、よく調べれば何とかなると思います。

Fortran ソース・プログラム

x64 RELEASE mode で DLL ライブラリにコンパイルする。

    module m_mandel
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, parameter :: kd = real64
      integer, parameter :: maxiter = 255
    contains
      subroutine mandel(nx, x, ny, y, m) bind(c, name = 'mandel')
        !DEC$ ATTRIBUTES DLLEXPORT :: mandel
        integer, value :: nx, ny
        real(kd), intent(in) :: x(nx), y(ny)
        integer, intent(out) :: m(ny, nx)
        integer :: ix, iy
        forall(ix = 1:nx, iy = 1:ny) m(iy, ix) = mandel_sub( cmplx(x(ix), y(iy), kind = kd) )
        return
      end subroutine mandel
    
      pure elemental integer function mandel_sub(c)
        complex(kd), intent(in) :: c
        complex(kd) :: z
        z = c
        do mandel_sub = maxiter, 1, -1
          if (abs(z) > 2.0_kd) exit
          z = z**2 + c
        end do
        return
      end function mandel_sub
    end module m_mandel
    

Python 側呼び出しプログラム

from ctypes import *
my_mod = CDLL("./x64/release/mandel3.dll")
my_mod.mandel.argtypes = [
   c_int32, np.ctypeslib.ndpointer(dtype=np.float64),
   c_int32, np.ctypeslib.ndpointer(dtype=np.float64),
   np.ctypeslib.ndpointer(dtype=np.int32)
   ]
my_mod.mandel.restype = c_void_p

x = linspace(-2.0, 2.0, 1024)
y = linspace(-2.0, 2.0, 1024)

m = np.zeros([y.size, x.size], int32)
my_mod.mandel(x.size, x, y.size, y, m)  


from PIL import Image
img = Image.new("RGB", (x.size, y.size))
for iy in range(y.size):
    for ix in range(x.size):
        k = m[iy, ix]
        img.putpixel((ix, iy), ((5 * k) % 256, k, (10 * k) % 256))


imshow(img)
img.save("mandel.png", "PNG")