Merge pull request #5 from gjbex/feature/f2py_omp

Feature/f2py omp

Author: gjbex

source-code/interfaciing-c-c++-fortran/F2py/Julia/Makefile

@@ -0,0 +1,11 @@
+F2PY = f2py3
+FFLAGS = -fopenmp -O2 -g
+VERSION = cpython-38-x86_64-linux-gnu
+
+all: julia.cpython-38-x86_64-linux-gnu.so
+
+julia.$(VERSION).so: julia_set_mod.f90 julia_set_omp_mod.f90
+	$(F2PY) -c -m julia --f90flags="$(FFLAGS)" $^
+
+clean:
+	$(RM) $(wildcard *.o) $(wildcard *mod) $(wildcard *.so)

source-code/interfaciing-c-c++-fortran/F2py/Julia/README.md

@@ -0,0 +1,15 @@
+# Julia set
+
+Illustration of using f2py to interface with OpenMP code.
+
+## What is it?
+
+1. `julia_set_mod.f90`: Fortran module defining a function that computes a
+   single point of a Julia set.
+
+1. `julia_set_omp_mod.f90`: Fortran module that computes a Julia set
+   parallelized with OpenMP.
+1. `julia_set.py`: Python script that initlalizes the arrays and computes
+   the Julia set using the Fortran implementation.
+1. `Makefile`: make file to build the shared library.
+1. `show_julia_set.py`: Python script to visualize the result.

source-code/interfaciing-c-c++-fortran/F2py/Julia/julia_set.py

@@ -0,0 +1,34 @@
+#!/usr/bin/env python
+
+from argparse import ArgumentParser
+import numpy as np
+import sys
+
+from julia import julia_set_omp_mod
+
+
+def init_z_values(nr_points):
+    real = np.linspace(-1.8, 1.8, nr_points)
+    imag = np.linspace(-1.8, 1.8, nr_points)
+    Real, Imag = np.meshgrid(real, imag)
+    return np.asfortranarray(Real + Imag*1.0j)
+
+
+def init_n_values(nr_points):
+    return np.zeros((nr_points, nr_points), dtype=np.int32, order='F')
+
+
+def main():
+    arg_parser = ArgumentParser(description='compute julia set')
+    arg_parser.add_argument('--n', type=int, default=100,
+                            help='number of points')
+    options = arg_parser.parse_args()
+    z = init_z_values(options.n)
+    n = init_n_values(options.n)
+    julia_set_omp_mod.julia_compute(z, n)
+    np.savetxt(sys.stdout, n, fmt='%4d')
+
+
+if __name__ == '__main__':
+    status = main()
+    sys.exit(status)

source-code/interfaciing-c-c++-fortran/F2py/Julia/julia_set_mod.f90

@@ -0,0 +1,26 @@
+module julia_set_mod
+    implicit none
+
+    private
+    integer, parameter :: max_iters = 255
+    real(kind=8), parameter :: max_norm = 2.0D00
+    complex(kind=8), parameter :: C = cmplx(-0.622772D00, 0.42193D00, kind=8)
+
+    public :: julia_iterate
+
+contains
+
+    elemental function julia_iterate(z0) result(n)
+        implicit none
+        complex(kind=8), intent(in) :: z0
+        integer :: n
+        complex(kind=8) :: z
+        z = z0
+        n = 0
+        do while (abs(z) < max_norm .and. n < max_iters)
+            z = z**2 + C
+            n = n + 1
+        end do
+    end function julia_iterate
+
+end module julia_set_mod

source-code/interfaciing-c-c++-fortran/F2py/Julia/julia_set_omp_mod.f90

@@ -0,0 +1,24 @@
+module julia_set_omp_mod
+    use julia_set_mod
+    implicit none
+
+    public :: julia_compute
+
+contains
+
+    subroutine julia_compute(z, n)
+        implicit none
+        complex(kind=8), dimension(:, :), intent(in) :: z
+        integer, dimension(:, :), intent(inout) :: n
+        integer :: i, j
+
+        !$omp parallel do default(none) private(i) shared(z, n) schedule(runtime)
+        do j = 1, size(z, 2)
+            do i = 1, size(z, 1)
+                n(i, j) = julia_iterate(z(i, j))
+            end do
+        end do
+        !$omp end parallel do
+    end subroutine julia_compute
+
+end module julia_set_omp_mod

source-code/interfaciing-c-c++-fortran/F2py/Julia/show_julia_set.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+from argparse import ArgumentParser
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+
+
+def read_data(file_name):
+    data = np.genfromtxt(file_name, dtype=np.uint8)
+    return data
+
+
+def visualize_data(data):
+    plt.imshow(np.log1p(data/255.0))
+    plt.axis('off')
+
+
+def main():
+    arg_parser = ArgumentParser(description='visualize julia set')
+    arg_parser.add_argument('file', help='file name')
+    options = arg_parser.parse_args()
+    data = read_data(options.file)
+    visualize_data(data)
+    plt.show()
+    return 0
+
+if __name__ == '__main__':
+    status = main()
+    sys.exit(status)

source-code/interfaciing-c-c++-fortran/F2py/README.md

@@ -14,3 +14,4 @@ What is it?
     is used to create a shared library and a wrapper.
 * `comp_pi_f2py.py`: test program to compute pi using a given number
     of iterations.
+* `Julia`: example of using Fortran OpenMP procedures from Python.