fortran90によるヒープソートとバブルソート

fortran90で配列のソートを行うプログラムを載せます。
クイックソート(ここでは載せません)はヒープソートよりも若干遅くなってしまいました。
プログラムの最適化が出来ていない、Fortranの性質(ルーチンの呼び出しが遅い)からかもしれません。

どういうソート方法なのかはwikipediaをご参照ください(ヒープソートバブルソート)。

計算時間はヒープソートの方が圧倒的に早いです。

クイックソートの方がヒープソートよりも早いです。コードは最速のクイックソートへ。

倍精度型の、サイズNの配列arrayを昇順(小さい数から大きい数へ)に並べ替えます。

バブルソート

平均計算量\(O(N^2)\)

subroutine bubblesort(N,array)
  !sikinote, 2016/08/08
  implicit none
  integer,intent(in)::N
  double precision,intent(inout)::array(1:N)
  integer::i,j
  double precision::t
   
  do i=1,N-1
     do j=i+1,N
        if(array(i) .gt. array(j))then
           t=array(i)
           array(i)=array(j)
           array(j)=t
        end if
     end do
  end do

  return
end subroutine bubblesort

ヒープソート

平均計算量\(O(N\log N)\)

subroutine heapsort(n,array)
  implicit none
  integer,intent(in) :: n
  double precision,intent(inout) :: array(1:n)
 
  integer ::i,k,j,l
  double precision :: t
 
  if(n.le.0)then
     write(6,*)"Error, at heapsort"; stop
  endif
  if(n.eq.1)return

  l=n/2+1
  k=n
  do while(k.ne.1)
     if(l.gt.1)then
        l=l-1
        t=array(L)
     else
        t=array(k)
        array(k)=array(1)
        k=k-1
        if(k.eq.1) then
           array(1)=t
           exit
        endif
     endif
     i=l
     j=l+l
     do while(j.le.k)
        if(j.lt.k)then
           if(array(j).lt.array(j+1))j=j+1
        endif
        if (t.lt.array(j))then
           array(i)=array(j)
           i=j
           j=j+j
        else
           j=k+1
        endif
     enddo
     array(i)=t
  enddo

  return
end subroutine heapsort

[adsense1]

配列arrayをどのように入れ替えたのか?という情報が欲しければ

バブルソート

平均計算量\(O(N^2)\)

subroutine bubblesort2(N,data,turn)
  !sikinote, 2016/08/08
  integer,intent(in)::N
  integer,intent(out)::turn(1:N)
  double precision,intent(inout)::data(1:N)
  integer::i,j,ti
  double precision::tmp

  do i=1,N
     turn(i)=i
  enddo

  do i=1,N-1
     do j=i+1,N
        if(data(i) .gt. data(j))then
           tmp=data(i)
           data(i)=data(j)
           data(j)=tmp
         
           ti=turn(i)
           turn(i)=turn(j)
           turn(j)=ti
        end if
     end do
  end do

  return
end subroutine bubblesort2

ヒープソート

平均計算量\(O(N\log N)\)

subroutine heapsort2(n,array,turn)
  implicit none
  integer,intent(in)::n
  integer,intent(out)::turn(1:n)
  double precision,intent(inout)::array(1:n)
 
  integer::i,k,j,l,m
  double precision::t
 
  if(n.le.0)then
     write(6,*)"Error, at heapsort"; stop
  endif
  if(n.eq.1)return

  do i=1,N
     turn(i)=i
  enddo

  l=n/2+1
  k=n
  do while(k.ne.1)
     if(l.gt.1)then
        l=l-1
        t=array(l)
        m=turn(l)
     else
        t=array(k)
        m=turn(k)
        array(k)=array(1)
        turn(k)=turn(1)
        k=k-1
        if(k.eq.1) then
           array(1)=t
           turn(1)=m
           exit
        endif
     endif
     i=l
     j=l+l
     do while(j.le.k)
        if(j.lt.k)then
           if(array(j).lt.array(j+1))j=j+1
        endif
        if (t.lt.array(j))then
           array(i)=array(j)
           turn(i)=turn(j)
           i=j
           j=j+j
        else
           j=k+1
        endif
     enddo
     array(i)=t
     turn(i)=m
  enddo

  return
end subroutine heapsort2

というプログラムで実装できます。

プログラム例1 ソートする


実際に実行してみます。
乱数で振った配列aを昇順に並べ替えるには

program main
  implicit none
  integer::i,N
  double precision,allocatable::a(:)
  double precision::t
 
  N=10
  allocate(a(0:N-1)); a=0d0

  do i=0,N-1
     call random_number(t)
     a(i)=t
     write(6,*)a(i)
  enddo

  call heapsort(size(a),a)
  do i=0,N-1
     write(6,*)i,a(i)
  enddo
 
  stop
end program main
   
subroutine heapsort(n,array)
  implicit none
  integer,intent(in) :: n
  double precision,intent(inout) :: array(1:n)
 
  integer ::i,k,j,l
  double precision :: t
 
  if(n.le.0)then
     write(6,*)"Error, at heapsort"; stop
  endif
  if(n.eq.1)return

  l=n/2+1
  k=n
  do while(k.ne.1)
     if(l.gt.1)then
        l=l-1
        t=array(L)
     else
        t=array(k)
        array(k)=array(1)
        k=k-1
        if(k.eq.1) then
           array(1)=t
           exit
        endif
     endif
     i=l
     j=l+l
     do while(j.le.k)
        if(j.lt.k)then
           if(array(j).lt.array(j+1))j=j+1
        endif
        if (t.lt.array(j))then
           array(i)=array(j)
           i=j
           j=j+j
        else
           j=k+1
        endif
     enddo
     array(i)=t
  enddo

  return
end subroutine heapsort

とすればよいです。実行結果は

$ gfortran main.f90
$ ./a.out
  0.99755959009261719    
  0.56682470761127335    
  0.96591537549612494    
  0.74792768547143218    
  0.36739089737475572    
  0.48063689875473148    
   7.3754263633984518E-002
   5.3552292777272470E-003
  0.34708128851801545    
  0.34224381607283505    
           0   5.3552292777272470E-003
           1   7.3754263633984518E-002
           2  0.34224381607283505    
           3  0.34708128851801545    
           4  0.36739089737475572    
           5  0.48063689875473148    
           6  0.56682470761127335    
           7  0.74792768547143218    
           8  0.96591537549612494    
           9  0.99755959009261719    
$

となります。

[adsense2]

プログラム例2 どのように並べ替えたかの順番を知る


program main
  implicit none
  integer::i,N
  integer,allocatable::turn(:)
  double precision,allocatable::a(:)
  double precision::t
 
  N=10
  allocate(a(0:N-1)); a=0d0
  allocate(turn(0:N-1)); turn=0

  do i=0,N-1
     call random_number(t)
     a(i)=t
     write(6,*)i,a(i)
  enddo
  write(6,*)

  call heapsort2(size(a),a,turn)
  do i=0,N-1
     write(6,*)turn(i),a(i)
  enddo
 
  stop
end program main
   
subroutine heapsort2(n,array,turn)
  implicit none
  integer,intent(in)::n
  integer,intent(out)::turn(1:n)
  double precision,intent(inout)::array(1:n)
 
  integer::i,k,j,l,m
  double precision::t
 
  if(n.le.0)then
     write(6,*)"Error, at heapsort"; stop
  endif
  if(n.eq.1)return

  do i=1,N
     turn(i)=i
  enddo

  l=n/2+1
  k=n
  do while(k.ne.1)
     if(l.gt.1)then
        l=l-1
        t=array(l)
        m=turn(l)
     else
        t=array(k)
        m=turn(k)
        array(k)=array(1)
        turn(k)=turn(1)
        k=k-1
        if(k.eq.1) then
           array(1)=t
           turn(1)=m
           exit
        endif
     endif
     i=l
     j=l+l
     do while(j.le.k)
        if(j.lt.k)then
           if(array(j).lt.array(j+1))j=j+1
        endif
        if (t.lt.array(j))then
           array(i)=array(j)
           turn(i)=turn(j)
           i=j
           j=j+j
        else
           j=k+1
        endif
     enddo
     array(i)=t
     turn(i)=m
  enddo

  return
end subroutine heapsort2

実行例

$ gfortran main.f90
$ ./a.out
           0  0.99755959009261719    
           1  0.56682470761127335    
           2  0.96591537549612494    
           3  0.74792768547143218    
           4  0.36739089737475572    
           5  0.48063689875473148    
           6   7.3754263633984518E-002
           7   5.3552292777272470E-003
           8  0.34708128851801545    
           9  0.34224381607283505    

           8   5.3552292777272470E-003
           7   7.3754263633984518E-002
          10  0.34224381607283505    
           9  0.34708128851801545    
           5  0.36739089737475572    
           6  0.48063689875473148    
           2  0.56682470761127335    
           4  0.74792768547143218    
           3  0.96591537549612494    
           1  0.99755959009261719

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です