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
!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
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
!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
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
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
$
$ ./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
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
$ ./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
