원주율 계산 (pi calculations) [병렬 계산] by 바죠

병렬 계산에서는 CPU 사용 갯수에 의존하지 않고 항상 같은 결과를 만들어야한다.
특히, random number를 사용하는 경우,  많은 사람들이 이를 무시하는 결과를 내어 놓는다.
하지만, 그것은 제대로 된 병렬 구현이 아니다.
랜덤 넘버 제너레이션 씨드가 주어지면 노드 수에 상관없이 단 하나의 결과를 주어야 한다.
다시말해서 원할 경우, 계산 결과를 정확하게 재생산 할 수 있어야 한다.
다른관점에서 보았을 때,
이것은 또한 랜덤 넘버가 노드별로 상관관계를 가질 수 있기도 하여 매우 주의해야 하는 항목이 된다.
http://en.wikipedia.org/wiki/Pi

정사각형에 랜덤 넘버들을 뿌려놓는다. 그 다음 원 안에 있는 것들의 비율을 계산함으로써 pi를 계산할 수 있다.
실제는 랜덤 넘버 두 개를 이용해서 원 안에 들어가는지 아닌지를 확인한다.
거의 완전히 병렬화가 가능한 알고리듬이다. 마지막 순간에만 노드들 사이의 통신이 필요한 수준이다.


4.d0*atan(1.d0)

그림 출처: google image

            2  processes are alive
    3.141701960000000    
    3.906250000000000    
            4  processes are alive
    3.141701960000000    
    1.957031250000000    
            8  processes are alive
    3.141701960000000    
    1.015625000000000    
            1  process is alive
    3.141701960000000    
    8.019531250000000    

cat input
100   nn
11 22 iseed1,iseed2


아래에 주어진 코드는 MPI 코드이다. 노드수에 의존하지 않은 결과를 주는 하나의 병렬 프로그램이다.
아래와 같이 노드수가 바뀌더라도 계산된 pi 값은 같다.
아울러 계산 시간도 함께 프린트해보았다.

seed가 주어지면 계산 결과는 반드시 결정되어야 한다. 물론,  병렬 계산에서는  "+ 순서"가  병렬 계산으로 인하여
바뀌기 때문에 이에 따른 아주 미세한 차이는 나올 수 있습니다. 예를 들어 대략 10^(-12) 정도 차이는 있을 수 있다.

아래와 같이 다소 차이가 나는 경우가 발생한다. 이 경우 완전히 같은 결과를 주어야 하나, 덧셈의 순서가 달라서 발생하는 오차라고 본다. 뒤에 나오는 심프슨 공식을 구간별로 적용한 경우에는 보다 노드 수에 따라서 보다 더 큰 편차를 주는 것을 알 수 있다.

 cat out11 out2 out4 out8
            1  process is alive
 number of intervals:
    3.141592653590426    
    1.070312500000000       sec
            2  processes are alive
 number of intervals:
    3.141592653589909    
   0.4296875000000000       sec
            4  processes are alive
 number of intervals:
    3.141592653589683    
   0.2148437500000000       sec
            8  processes are alive
 number of intervals:
    3.141592653589815    
   0.1367187500000000       sec


------------------------------------------------------------------


!234567890
!      Written by In-Ho Lee, KRISS, February, 2011
       program pi_parallel
       implicit none
       include 'mpif.h'
       integer ierr,kount,iroot
       integer i,jdum,kdum,nn,iseed1,iseed2,ij,kl
       real*8 tmp,tmr,r,x,y, t1,t2
       real*8 tmp0,tmr0
       real ranmar
       integer myid,nproc
       integer istart,ifinish,n1,n2

       call MPI_INIT( ierr )
       call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
       call MPI_COMM_SIZE( MPI_COMM_WORLD, nproc, ierr )
       if(myid == 0 .and. nproc > 1) print *,  nproc," processes are alive"
       if(myid == 0 .and. nproc ==1) print *,  nproc," process is alive"

       t1=MPI_Wtime()

       if(myid == 0) then     ! -------=== { process id =0
       iseed1=11
       iseed2=22
       nn=100

       read(5,*) nn
       read(5,*) iseed1,iseed2
                     endif    ! -------=== } process id =0
       iroot=0 ; kount=1
       call MPI_BCAST(nn, kount,MPI_INTEGER,iroot,MPI_COMM_WORLD,ierr)
       call MPI_BCAST(iseed1, kount,MPI_INTEGER,iroot,MPI_COMM_WORLD,ierr)
       call MPI_BCAST(iseed2, kount,MPI_INTEGER,iroot,MPI_COMM_WORLD,ierr)

       tmp=0.d0
       tmr=0.d0

       n1=1
       n2=nn
       call equal_load(n1,n2,nproc,myid,istart,ifinish)
!      print*, istart,ifinish
!      do i=1,nn
       do i=istart,ifinish

       ij=iseed1 +i      ; kl=iseed2
       ij=mod(ij, 31328 +1) ; kl=mod(kl, 30081 +1)
!      print*, ij,kl
       call rmarin(ij,kl)

       do jdum=1,1000
       do kdum=1,1000

       x=ranmar()-0.5
       y=ranmar()-0.5
       r=sqrt(x*x+y*y)
       tmr=tmr+1.d0
       if(r < 0.5d0) then
       tmp=tmp+1.d0
                     endif

       enddo
       enddo
       enddo
       t2=MPI_Wtime()
       iroot=0 ; kount=1
       call MPI_REDUCE(tmp,tmp0,kount,MPI_DOUBLE_PRECISION,MPI_SUM,iroot,MPI_COMM_WORLD,ierr)
       call MPI_REDUCE(tmr,tmr0,kount,MPI_DOUBLE_PRECISION,MPI_SUM,iroot,MPI_COMM_WORLD,ierr)


       if(myid == 0) then     ! -------=== { process id =0
       print*, tmp0/tmr0*4.d0
       print*, t2-t1
                     endif    ! -------=== } process id =0

       call MPI_FINALIZE(ierr)
       stop
       end program pi_parallel



      subroutine rmarin(ij,kl)
!  This subroutine and the next function generate random numbers. See
!  the comments for SA for more information. The only changes from the
!  orginal code is that (1) the test to make sure that RMARIN runs first
!  was taken out since SA assures that this is done (this test didn't
!  compile under IBM's VS Fortran) and (2) typing ivec as integer was
!  taken out since ivec isn't used. With these exceptions, all following
!  lines are original.

! This is the initialization routine for the random number generator
!     RANMAR()
! NOTE: The seed variables can have values between:    0 <= IJ <= 31328
!                                                      0 <= KL <= 30081
      real u(97), c, cd, cm
      integer i97, j97
      common /raset1/ u, c, cd, cm, i97, j97
      if( ij .lt. 0  .or.  ij .gt. 31328  .or.  &
          kl .lt. 0  .or.  kl .gt. 30081 ) then
          print '(a)', ' The first random number seed must have a value  &
    & between 0 and 31328'
          print '(a)',' The second seed must have a value between 0 and  &
    & 30081'
            stop
      endif
      i = mod(ij/177, 177) + 2
      j = mod(ij    , 177) + 2
      k = mod(kl/169, 178) + 1
      l = mod(kl,     169)
      do 2 ii = 1, 97
         s = 0.0
         t = 0.5
         do 3 jj = 1, 24
            m = mod(mod(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = mod(53*l+1, 169)
            if (mod(l*m, 64) .ge. 32) then
               s = s + t
            endif
            t = 0.5 * t
 3       continue
         u(ii) = s
 2    continue
      c = 362436.0 / 16777216.0
      cd = 7654321.0 / 16777216.0
      cm = 16777213.0 /16777216.0
      i97 = 97
      j97 = 33
      return
      end
      function ranmar()
      real u(97), c, cd, cm
      integer i97, j97
      common /raset1/ u, c, cd, cm, i97, j97
         uni = u(i97) - u(j97)
         if( uni .lt. 0.0 ) uni = uni + 1.0
         u(i97) = uni
         i97 = i97 - 1
         if(i97 .eq. 0) i97 = 97
         j97 = j97 - 1
         if(j97 .eq. 0) j97 = 97
         c = c - cd
         if( c .lt. 0.0 ) c = c + cm
         uni = uni - c
         if( uni .lt. 0.0 ) uni = uni + 1.0
         ranmar = uni
      return
      end
       subroutine equal_load(n1,n2,nproc,myid,istart,ifinish)
!      Written by In-Ho Lee, KRISS, September (2006)
       implicit none
       integer nproc,myid,istart,ifinish,n1,n2
       integer iw1,iw2
       iw1=(n2-n1+1)/nproc ; iw2=mod(n2-n1+1,nproc)
       istart=myid*iw1+n1+min(myid,iw2)
       ifinish=istart+iw1-1 ; if(iw2 > myid) ifinish=ifinish+1
!      print*, n1,n2,myid,nproc,istart,ifinish
       if(n2 < istart) ifinish=istart-1
       return
       end


그림 출처: google image


 

------------------------------------------------------------------

PS.



!      Written by In-Ho Lee, KRISS, February (2011)
       program pi_parallel2
       implicit none
       include 'mpif.h'
       real*8  f, a, dx, x, asum, tmp, pi
       integer nx, ix
       integer n1,n2,istart,ifinish
       integer ierr, nproc, myid
       real*8 t1,t2


       f(a) = 4.d0 / (1.d0 + a*a)


       call MPI_INIT(ierr)
       call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ierr)
       call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
       if(myid == 0 .and. nproc > 1) print *,  nproc," processes are alive"
       if(myid == 0 .and. nproc ==1) print *,  nproc," process is alive"

       if(myid == 0) then     ! -------=== { process id =0
       print*, 'number of intervals:'
       read(*,*) nx
                     endif    ! -------=== } process id =0

       t1=MPI_Wtime()
       call MPI_BCAST(nx,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       dx = 1.0d0 / dfloat(nx)

       asum = 0.0d0
       n1=1 ; n2=nx
       call equal_load(n1,n2,nproc,myid,istart,ifinish)
       do 10 ix = istart, ifinish
         x = dx*(dfloat(ix)-0.5d0)
         asum = asum + f(x)
  10   continue
       tmp = dx*asum

       call MPI_REDUCE(tmp,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
       if(myid == 0) then     ! -------=== { process id =0
       print*, pi
                     endif    ! -------=== } process id =0

       t2=MPI_Wtime()
       if(myid == 0) then     ! -------=== { process id =0
       print*, t2-t1,' sec '
                     endif    ! -------=== } process id =0

       call MPI_FINALIZE(ierr)
       stop
       end program pi_parallel2
!234567890
       subroutine equal_load(n1,n2,nproc,myid,istart,ifinish)
!      Written by In-Ho Lee, KRISS, September (2006)
       implicit none
       integer nproc,myid,istart,ifinish,n1,n2
       integer iw1,iw2
       iw1=(n2-n1+1)/nproc ; iw2=mod(n2-n1+1,nproc)
       istart=myid*iw1+n1+min(myid,iw2)
       ifinish=istart+iw1-1 ; if(iw2 > myid) ifinish=ifinish+1
!      print*, n1,n2,myid,nproc,istart,ifinish
       if(n2 < istart) ifinish=istart-1
       return
       end

모든 계산 노드에서 위의 서브루틴을 부르면, 각 노드에서 계산의 영역을 반환해 준다.
반환되는 값은 istart,ifinish이다.
인풋은 myid,nproc, n1,n2이다.
즉, n1,n2는 노드의 수가 1개일 때, 계산하는 영역이다.
이것을 nproc의 노드가 나누어서 일을 수행할 수 있도록 do loop의 시작과 끝을 계산 노드별로 계산을 해 준다. 노드의 숫자에 땨라서 균등하게 나누어지지 않는 경우가 발생할 수 있다.  하지만, 전체 계산의 영역을 하나도 빠짐없이 수행할 수 있게 해준다. 즉, 노드 수가 1개 일때와 차이가 없이 완전히  전 영역을 커버하게 해 준다.


정밀도 부문에서 두 번째 프로그램이 더 뛰어난것 같다.

a.out
            1  process is alive
 number of intervals:
10000000
    3.141592653589731    
   0.1054687500000000       sec
FORTRAN STOP


이번에는 좀 더 체계적인 방법으로 정적분을 수행해 본다.
이 경우, 아래와 같이 심슨의 룰을 활용한다.  아래의 구현은 좀 더 일반적인 1차원 정적분에 응용할 수 있는 코드에 해당한다.
함수 func을 재 정의함으로써 보다 더 일반적으로 프로그램을 활용할 수 있게 만들었다.
Simpson rule을 사용하는 방법:
http://incredible.egloos.com/4494526
http://en.wikipedia.org/wiki/Simpson_rule

------------------------------------------------------------------


!      Written by In-Ho Lee, KRISS, February (2011)
       program pi_parallel3
       implicit none
       include 'mpif.h'
       real*8  f, x, tmp, pi
       integer nx, ix
       integer nn1,nn2,istart,ifinish,m1
       integer ierr, nproc, myid
       real*8 t1,t2
       real*8 aa,bb
       real*8 a1,b1
       external func


       call MPI_INIT(ierr)
       call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ierr)
       call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
       if(myid == 0 .and. nproc > 1) print *,  nproc," processes are alive"
       if(myid == 0 .and. nproc ==1) print *,  nproc," process is alive"

       if(myid == 0) then     ! -------=== { process id =0
       print*, 'Simpson s rule'
       print*, 'number of intervals:'
       read(*,*) nx
                     endif    ! -------=== } process id =0

       t1=MPI_Wtime()
       call MPI_BCAST(nx,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       aa=0.d0
       bb=1.d0

       nn1=1 ; nn2=nx
       call equal_load(nn1,nn2,nproc,myid,istart,ifinish)

       m1=0
       do ix=istart,ifinish
       m1=m1+1
       enddo
       a1=aa+(bb-aa)/float(nx-1)*float(istart-1)
       b1=aa+(bb-aa)/float(nx-1)*float(ifinish-1)
       call simpson(func,m1,a1,b1,tmp)

       call MPI_REDUCE(tmp,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
       if(myid == 0) then     ! -------=== { process id =0
       print*, pi
                     endif    ! -------=== } process id =0

       t2=MPI_Wtime()
       if(myid == 0) then     ! -------=== { process id =0
       print*, t2-t1,' sec '
                     endif    ! -------=== } process id =0

       call MPI_FINALIZE(ierr)
       stop
       end program pi_parallel3


!234567890
       subroutine equal_load(n1,n2,nproc,myid,istart,ifinish)
!      Written by In-Ho Lee, KRISS, September (2006)
       implicit none
       integer nproc,myid,istart,ifinish,n1,n2
       integer iw1,iw2
       iw1=(n2-n1+1)/nproc ; iw2=mod(n2-n1+1,nproc)
       istart=myid*iw1+n1+min(myid,iw2)
       ifinish=istart+iw1-1 ; if(iw2 > myid) ifinish=ifinish+1
!      print*, n1,n2,myid,nproc,istart,ifinish
       if(n2 < istart) ifinish=istart-1
       return
       end
       subroutine simpson(f,n,aa,bb,rslt)
!      Written by In-Ho Lee, KRISS, 2006
       implicit none
       real*8 f
       integer n
       real*8 rslt,aa,bb
       real*8 h,xx
       integer j
       logical lodd
       integer m

       m=2*n
       h=(bb-aa)/float(m)
       rslt=(f(aa)+f(bb))
       lodd=.true.
       do j=1,m-1
        xx=aa+h*float(j)
       if(lodd)then
       rslt=rslt+4.0d0*f(xx)
       else
         rslt=rslt+2.0d0*f(xx)
       endif
       lodd=(.not. lodd)
       enddo

       rslt=rslt*h/3.0d0
       return
       end

       real*8 function func(x)
       implicit none
       real*8 x

       func = 4.d0 / (1.d0 + x*x)
       return
       end

a.out<input
            1  process is alive
 Simpson s rule
 number of intervals:
    3.141592653588563    
    6.316406250000000       sec
FORTRAN STOP

4  processes are alive
 Simpson s rule
 number of intervals:
    3.141592653589751    
    1.593750000000000       sec

http://en.wikipedia.org/wiki/Definite_integral
http://en.wikipedia.org/wiki/Simpson_rule



grep sec out11 out2 out4 out8
out11:    6.316406250000000       sec
out2:    3.171875000000000       sec
out4:    1.589843750000000       sec
out8:   0.7929687500000000       sec
[ihlee@helix pi_mpi]$ bc -l
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'.
6.316/0.792
7.97474747474747474747
6.316/1.58
3.99746835443037974683
6.316/3.1718
1.99129831641339302604

MPI_Wtime()을 사용해서 wall-clock 시간을 초단위로 잴 수 있다. 이것은 병렬 효율성을 체크할 때 필요하다.
지금 풀고 있는 문제가 소위
embarrassingly parallel problem 이라는 것을 보여 주는 그래프를 아래에 보였다. 노드 수에 거의 정비례하는 speedup. 
http://en.wikipedia.org/wiki/Embarrassingly_parallel

cat out11 out2 out4 out8
1 process is alive
Simpson s rule
number of intervals:
3.141592653588563
6.316406250000000 sec
2 processes are alive
Simpson s rule
number of intervals:
3.141592653590160
3.171875000000000 sec
4 processes are alive
Simpson s rule
number of intervals:
3.141592653589751
1.589843750000000 sec
8 processes are alive
Simpson s rule
number of intervals:
3.141592614205191
0.7929687500000000 sec

-----------------------------------------


새로운 랜덤 넘버 제너레이터(포트란 90):
http://incredible.egloos.com/4494190
아래 모듈 출처:
http://www.sdsc.edu/~tkaiser/f90.html
module ran_mod
! ran1 returns a uniform random number between 0-1
! the seed is optional and used to reset the generator
contains
function ran1(my_seed)
implicit none
real(8) ran1,r
integer, optional ,intent(in) :: my_seed ! optional argument not changed in the routine
integer,allocatable :: seed(:)
integer the_size,j
if(present(my_seed))then ! use the seed if present
call random_seed(size=the_size) ! how big is the intrisic seed?
allocate(seed(the_size)) ! allocate space for seed
do j=1,the_size ! create the seed
seed(j)=abs(my_seed)+(j-1) ! abs is generic
enddo
call random_seed(put=seed) ! assign the seed
deallocate(seed) ! deallocate space
endif
call random_number(r)
ran1=r
end function ran1
end module


program darwin
use ran_mod ! interface required if we have
! optional or intent arguments
real(8) x,y
x=ran1(my_seed=12345) ! we can specify the name of the argument
y=ran1()
write(*,*) x,y
x=ran1(12345) ! with only one optional argument we don't need to
y=ran1()
write(*,*) x,y
end program


a.out
2.9463771517441728E-003 2.9468539889592194E-003
2.9463771517441728E-003 2.9468539889592194E-003


http://incredible.egloos.com/3021478

Multiple Precision Computation
다중 정밀도 계산을 수행하지 않는 한 배정밀도 보다 더 정밀한 계산을 할 수 없다.

씨리얼 코드의 구간별 시간 측정이 코드 최적화의 시작점이다.
http://incredible.egloos.com/3047318


---------------------------------------------------------
PS.

동일한 씨드의 중복 사용을 피하기 위해서 아래와 같이 점프를 시도할 수 있다. 
단순한 공식으로 랜덤 씨드를 만들어 낼 경우, 병렬 프로그램에서 다음 순서의 싸이클에서 뜻밖에 노드들에게 동일한 씨드가 들어갈 수 있다. 사용하는 데이터가 거의 수렴한 경우에는 이 경우 거의 같은 계산을 다시 수행하는 사태가 발생할 수 있다.  수렴했다고 착각할 수 있다.
 
이러한 경우를 피하기 위해서 0 번 노드에서 씨드 하나를 만들고 (중복을 피할 수 있게 굉장히 큰 스케일에서 하나 만든다. 씨드로 부터 결정론적으로 만들어진다.)  이것을 모든 노드에 분배하고, 각 노드에서는 이를 다시 이용해서 다시 노드별로, 알고리듬에 의해서, 새로운 씨드를 만들어내어서 사용한다.
 


       do jtry=1,mtry
       if(myid == 0) kdum=int(ran1(ksave+iseed+jtry)*1.d8)+1
       iroot=0 ; kount=1
       call MPI_BCAST(kdum, kount,MPI_INTEGER,iroot,MPI_COMM_WORLD,ierr)
      
       n1=1 ; n2=1024
       call equal_load(n1,n2,nproc,myid,istart,ifinish)
       do itry=istart,ifinish
       kdum=int(ran1(kdum+itry)*1.d8)+1 ; tmr=ran1(kdum)




-------------------------

포트란 90에서 지원하는 랜덤 넘버 사용하기:
우선 씨드를 주고 초기화한다. 그다음 사용한다.
nseed=222
call random_seed(nseed)
call random_number(tmp)

핑백

덧글

  • 필군 2011/02/24 01:29 # 답글

    저는 통상 random_seed()를 써서 초기화 했는데 다른 방식을 쓰시는 듯 하군요 - 프로세서간에 다른 initial seed를 주는게 관건인듯 합니다.
  • 바죠 2011/02/24 08:48 #

    덧글 감사합니다. 그렇습니다. 랜덤 넘버 제너레이터가 좀 특이한 케이스입니다.
  • 바죠 2011/02/24 13:48 #


    아래의 URL에서 복사하고 약간 수정한 루틴이 유용할 것 같습니다.

    random number module from:
    http://www.sdsc.edu/~tkaiser/f90.html

    !234567890
    module ran_mod
    ! ran1 returns a uniform random number between 0-1
    ! the seed is optional and used to reset the generator
    contains
    function ran1(my_seed)
    implicit none
    real(8) ran1,r
    integer, optional ,intent(in) :: my_seed ! optional argument not changed in the routine
    integer,allocatable :: seed(:)
    integer the_size,j
    if(present(my_seed))then ! use the seed if present
    call random_seed(size=the_size) ! how big is the intrisic seed?
    allocate(seed(the_size)) ! allocate space for seed
    do j=1,the_size ! create the seed
    seed(j)=abs(my_seed)+(j-1) ! abs is generic
    enddo
    call random_seed(put=seed) ! assign the seed
    deallocate(seed) ! deallocate space
    endif
    call random_number(r)
    ran1=r
    end function ran1
    end module

  • 2011/02/24 13:26 # 답글 비공개

    비공개 덧글입니다.
  • 2011/02/24 13:28 # 답글 비공개

    비공개 덧글입니다.
  • 2011/02/24 13:33 # 답글 비공개

    비공개 덧글입니다.
  • 바죠 2013/02/18 17:23 # 답글

    implicit none
    integer nseed
    real*8 tmp,tmq,aa
    real*8 histo(0:100000000)
    real*8 avg,sig,tmr
    integer i

    nseed=222
    call random_seed(nseed)
    call random_number(tmp)

    histo=0.d0

    avg=0.d0
    sig=0.d0
    tmr=0.d0

    aa=1.d0
    aa=100.d0
    do i=1,1000000
    call random_number(tmp)
    tmq=-aa*log(tmp)
    avg=avg+tmq
    sig=sig+tmq**2
    tmr=tmr+1.d0
    histo(int(tmq*10))=histo(int(tmq*10))+1.d0
    enddo

    avg=avg/tmr
    sig=sig/tmr
    sig=sqrt(sig-avg**2)
    print*, avg,sig

    open(9,file='fort.9',form='formatted')
    do i=0,10000
    write(9,*) i, histo(i)
    enddo
    close(9)

    stop
    end
댓글 입력 영역

최근 포토로그