[F2PY] 03. Advanced Example: Using OpenMP

Screen Shot 2018-04-23 at 1.50.11 PM.png

F2Py란?

수치 계산의 '우사인 볼트'라 할 수 있는 Fortran은 단점으로 다양한 파일의 입출력 같은 부가기능이 약합니다.
다재다능 Python은 단점으로 단순 계산이 조금 느립니다.
그래서 Python 뼈대에 계산 부분만 Fortran을 불러올 수 있게 도와주는 것이 F2PY입니다.
F2PY는 Numpy에 기본으로 내재되어 있습니다.
이 포스팅은 기본적으로 Linux 환경에서 Python3를 이용함을 알려드립니다.

1. Intro

OpenMP는 병렬화 Parallelizing를 시키는 가장 쉬운 방법입니다. OpenMP는 공유 메모리 시스템 Shared Memory System에 사용하기에 적합합니다. 최근 씨피유 CPU는 기본적으로 멀티쓰레드 Multi-Thread가 가능하기 때문에, 한 대의 개인용 컴퓨터 혹은 워크스테이션에서 멀티쓰레드를 이용하여 병렬화시킬 때 OpenMP가 유용합니다.

OpenMP의 공식 홈페이지는 여기입니다: http://www.openmp.org/
이외에도 각자의 언어에 맞는 OpenMP 사용법은 검색엔진에서 다수를 찾을 수 있습니다.

아래는 Fortran에서 멀티쓰레드를 생성하고 각각의 쓰레드 번호를 출력하는 기본 프로그램입니다.

module OTmod
  !$ use omp_lib
  implicit none

  public :: get_threads, set_num_threads
contains

  function get_threads() result(nt)
    integer :: myid,nthreads
    integer :: nt
    nt = 0
    !$ nt = omp_get_max_threads()

    !$OMP PARALLEL default(none) private(myid) &
    !$OMP shared(nthreads)

    ! Determine the number of threads and their id
    !$ myid = OMP_GET_THREAD_NUM()
    !$ nthreads = OMP_GET_NUM_THREADS()

    print*,nthreads,myid
    !$OMP END PARALLEL
  end function get_threads

  function set_num_threads(nn) result(nt)
    integer :: nn,nt
    !F2PY INTENT(IN) :: nn

    !$ nt = omp_get_max_threads()
    print*,"Current Max # of Threads=",nt
    !$ call omp_set_num_threads(nn)
    !$ nt = omp_get_max_threads()
    print*,"Current Max # of Threads=",nt
  end function set_num_threads

end module OTmod

위는 Fortran90 코드이며, StackOverflow의 한 페이지의 도움을 받았습니다.
https://stackoverflow.com/questions/32746104/difficulty-getting-openmp-to-work-with-f2py?rq=1

!$로 시작하는 부분이 OpenMP에서 해석하는 명령어 줄입니다.
("!"는 Fortran90에서 코멘트를 뜻합니다.)

OpenMP를 이용하기 위해서는 gfortran의 경우 -fopenmp 라는 옵션이 컴파일할 때 필요하며, f2py는 따로 -lgomp라는 옵션을 필요로 합니다. 따라서 위 프로그램의 컴파일 명령어는 다음과 같습니다.

$ f2py3 -m OTfor --fcompiler=gfortran --f90flags='-fopenmp' -lgomp -c OTmod.f90

그러면 "OTfor.cpython-34m.so"라는 파일이 생성되며, Python에서 읽을 준비가 되었습니다.

위 함수를 불러들이고 실행할 모체 Python 프로그램은 다음과 같습니다.

from OTfor import otmod

otmod.set_num_threads(5)
num=otmod.get_threads()
print("# of Threads = {}".format(num))

위에서 쓰레드의 갯수를 5개로 제한하였습니다.

위 프로그램을 실행시키면 다음과 같은 결과가 나옵니다.

 Current Max # of Threads=          16
 Current Max # of Threads=           5
           5           0
           5           4
           5           2
           5           1
           5           3
# of Threads = 5

5개의 쓰레드의 ID가 무작위 Random 순서로 찍혔습니다.

2. Speed by Number of Threads

멀티쓰레드가 시간을 얼마나 단축시키는지 알아보기 위하여 공간 상의 두 점 사이의 거리 Euclidean Distance를 계산하는 프로그램을 만들었습니다.

module calc_dist_omp_module
  !$ use omp_lib
  implicit none

  public :: calc_dist_2D,set_num_threads
contains

  function set_num_threads(nn) result(nt)
    integer :: nn,nt
    !F2PY INTENT(IN) :: nn
    !F2PY INTENT(OUT) :: nt
    
    !$ nt = omp_get_max_threads()
    print*,"Max # of Threads=",nt
    !$ call omp_set_num_threads(nn)
    !$ nt = omp_get_max_threads()
    print*, "Set Max # of Threads=",nt
  end function set_num_threads
 
  real(8) function calc_dist(aa,bb,nn) 
    integer :: nn,ii
    real(8),intent(in) :: aa(nn),bb(nn) 
    real(8) :: calc_sqdist
    !F2PY INTENT(HIDE) :: nn
    !F2PY depend(nn) aa
    !F2PY INTENT(OUT) :: calc_dist

    calc_sqdist=0.d0
    do ii=1,nn
       calc_sqdist=calc_sqdist+(aa(ii)-bb(ii))**2
    enddo
    calc_dist=sqrt(calc_sqdist)
    return
  end function calc_dist

  subroutine calc_dist_2D(A,B,dist1d,m,n)
  !!!
  !!! Get 2D array A(m,n) and B(m,n), and calculate distance along the 1st axis
  !!! Produce 1D array B(n) the value of which is distance A and B (m dimension)
  !!!
    integer :: m,n,i,j
    real(8), dimension(m,n),intent(in) :: A,B
    real(8), dimension(n) ,intent(out) :: dist1d
    integer :: myid,nthreads
!    real(8) :: calc_dist
    !F2PY INTENT(HIDE) :: m,n
    !F2PY depend(m,n) A
    !F2PY INTENT(OUT) :: dist1d

  !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(A,B,dist1d,m,n)

!    !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(A,B,dist1d,m,n)
!    ! Check T_id
!    !$ myid = OMP_GET_THREAD_NUM()
!    !$ nthreads = OMP_GET_NUM_THREADS()

!    print*,nthreads,myid
!    !$OMP DO
    do j=1,n
       dist1d(j)=calc_dist(A(:,j),B(:,j),m)
    enddo
!    !$OMP END DO
!    !$OMP END PARALLEL
  !$OMP END PARALLEL DO
  end subroutine calc_dist_2D

end module calc_dist_omp_module

위 프로그램의 컴파일은 다음과 같이 하고

$ f2py3 -c --fcompiler=gfortran --f90flags='-fopenmp' -lgomp -m calc_dist_omp calc_dist_mod.f90

위 함수를 불러들일 모체 Python 프로그램은 다음과 같습니다.

from calc_dist_omp import calc_dist_omp_module as cdistmod
import numpy as np
import timeit
import sys

#print("*** Module Document ***")
#print(cdistmod.__doc__)
#print("*** Function Document ***")
#print(cdistmod.calc_dist_2d.__doc__)

### Prepare input array for test
mm=100; nn=1000000
xx=np.linspace(0.,mm,num=nn,endpoint=False)
xx=np.tile(xx,mm).reshape([mm,nn])
array1=np.sin(xx)
array2=np.cos(xx)
print(array1.shape)

for max_threads in range(1,7,1):
  print("\n")
  cdistmod.set_num_threads(max_threads)
  print("*** Start Loop for {} ***".format(max_threads))
  time_record=[]
  for i in range(20):
    start = timeit.default_timer()
    dist=cdistmod.calc_dist_2d(array1,array2)
    stop = timeit.default_timer()
    time_record.append(stop-start)
  print("1st: {:.5f}sec 2nd: {:.5f}sec 3rd: {:.5f}sec".format(*sorted(time_record)[:3]))
  print(dist.shape,dist[:5])

### Mimic with Numpy
array1=array1.T
array2=array2.T
print("\n*** Start Loop ***")
time_record=[]
for i in range(10):
  start = timeit.default_timer()
  dist=np.sqrt(np.power(array1-array2,2).sum(axis=1))
  stop = timeit.default_timer()
  time_record.append(stop-start)
print("1st: {:.5f}sec 2nd: {:.5f}sec 3rd: {:.5f}sec".format(*sorted(time_record)[:3]))
print(dist.shape,dist[:5])

위 프로그램은 [100, 1000000] 크기의 행렬 A와 B를 만들고, 두 행렬간의 첫번째 축 1st Axis를 따라 거리를 계산하여 [1,1000000] 크기의 거리 행렬 (혹은 벡터)를 생성합니다. 쓰레드 갯수를 1부터 6까지 변화시켜 같은 계산을 하는데 걸리는 시간 변화를 알아봅니다. 또한 끝에는 Numpy만을 이용하여 같은 결과를 생성할 때 걸리는 시간도 출력합니다. 위 프로그램을 실행시키면 다음과 같은 결과가 나옵니다.

(100, 1000000)


 Max # of Threads=          16
 Set Max # of Threads=           1
*** Start Loop for 1 ***
1st: 1.39093sec 2nd: 1.39195sec 3rd: 1.39217sec
(1000000,) [ 10.           9.99899995   9.9979998    9.99699955   9.9959992 ]


 Max # of Threads=           1
 Set Max # of Threads=           2
*** Start Loop for 2 ***
1st: 1.30177sec 2nd: 1.30178sec 3rd: 1.30185sec
(1000000,) [ 10.           9.99899995   9.9979998    9.99699955   9.9959992 ]


 Max # of Threads=           2
 Set Max # of Threads=           3
*** Start Loop for 3 ***
1st: 1.26300sec 2nd: 1.26310sec 3rd: 1.26315sec
(1000000,) [ 10.           9.99899995   9.9979998    9.99699955   9.9959992 ]


 Max # of Threads=           3
 Set Max # of Threads=           4
*** Start Loop for 4 ***
1st: 1.23960sec 2nd: 1.24023sec 3rd: 1.24043sec
(1000000,) [ 10.           9.99899995   9.9979998    9.99699955   9.9959992 ]


 Max # of Threads=           4
 Set Max # of Threads=           5
*** Start Loop for 5 ***
1st: 1.23640sec 2nd: 1.23802sec 3rd: 1.23905sec
(1000000,) [ 10.           9.99899995   9.9979998    9.99699955   9.9959992 ]


 Max # of Threads=           5
 Set Max # of Threads=           6
*** Start Loop for 6 ***
1st: 1.23485sec 2nd: 1.23492sec 3rd: 1.23599sec
(1000000,) [ 10.           9.99899995   9.9979998    9.99699955   9.9959992 ]

*** Start Loop ***
1st: 3.28598sec 2nd: 3.28677sec 3rd: 3.28678sec
(1000000,) [ 10.           9.99899995   9.9979998    9.99699955   9.9959992 ]

3. Summary and Conclusion

위에 나온 결과를 Google Sheet을 이용하여 간단한 그래프로 나타내보았습니다.
chart (1).png

길이 1,000,000인 2개의 Array 사이에 100차원의 거리를 구하는 정도의 계산이라면 굳이 쓰레드를 4보다 많이 이용할 필요는 없어보입니다. Loop를 나눠 일할 일꾼이 늘어나도 단축되는 시간에는 한계가 있습니다. 일을 나누고 각자 한 일을 합산하는 작업에 일정 부분 시간이 걸리기 때문입니다. 또한 병렬화가 불가능한 부분에서 생기는 병목현상도 고려해야 합니다.

OpenMP는 손쉽게 멀티쓰레드를 이용할 수 있도록 도와줍니다. 적당한 갯수의 쓰레드로 최적의 효율을 찾을 수 있습니다.

"""
제 개인적 목표는