RとFortranの演算速度の比較

そろそろFortranもお勉強しようかとおもい、プログラムをしてみることに!!

作業は、2地点間の距離と角度を求めるというシンプルなもの
ただ、データ数が多く、組み合わせが全部で22,232,390もある。

まずRでプログラミング
特に変なことはしていないと思いますが...

Tree<-read.csv("TreePlot",header=T)
IncW3<-read.table("IncW3Th10S3",header=T,sep=" ")
IncW5<-read.table("IncW5Th10S3",header=T,sep=" ")
IncW11<-read.table("IncW11Th10S3",header=T,sep=" ")

Tree<-as.matrix(Tree[,1:5])
IncW3<-as.matrix(IncW3)
IncW5<-as.matrix(IncW5)
IncW11<-as.matrix(IncW11)

DistW3<-matrix(rep(NA,nrow(Tree)*nrow(IncW3)*7),ncol=7)
DistW5<-matrix(rep(NA,nrow(Tree)*nrow(IncW5)*7),ncol=7)
DistW11<-matrix(rep(NA,nrow(Tree)*nrow(IncW11)*7),ncol=7)

for (i in 1:nrow(Tree)){
 for (j in 1:nrow(IncW3)){
 DistW3[i*j,1:4]<-c(Tree[i,4:5],IncW3[j,1:2])
 DistW3[i*j,5]<-sqrt((Tree[i,4]-IncW3[j,1])^2+(Tree[i,5]-IncW3[j,2])^2)
 DistW3[i*j,6]<-atan((Tree[i,5]-IncW3[j,2])/(Tree[i,4]-IncW3[j,1]))
 DistW3[i*j,7]<-atan((Tree[i,5]-IncW3[j,2])/(Tree[i,4]-IncW3[j,1]))*180/pi
 }
 for (j in 1:nrow(IncW5)){
 DistW5[i*j,1:4]<-c(Tree[i,4:5],IncW5[j,1:2])
 DistW5[i*j,5]<-sqrt((Tree[i,4]-IncW5[j,1])^2+(Tree[i,5]-IncW5[j,2])^2)
 DistW5[i*j,6]<-atan((Tree[i,5]-IncW5[j,2])/(Tree[i,4]-IncW5[j,1]))
 DistW5[i*j,7]<-atan((Tree[i,5]-IncW5[j,2])/(Tree[i,4]-IncW5[j,1]))*180/pi
 }
 for (j in 1:nrow(IncW11)){
 DistW11[i*j,1:4]<-c(Tree[i,4:5],IncW11[j,1:2])
 DistW11[i*j,5]<-sqrt((Tree[i,4]-IncW11[j,1])^2+(Tree[i,5]-IncW11[j,2])^2)
 DistW11[i*j,6]<-atan((Tree[i,5]-IncW11[j,2])/(Tree[i,4]-IncW11[j,1]))
 DistW11[i*j,7]<-atan((Tree[i,5]-IncW11[j,2])/(Tree[i,4]-IncW11[j,1]))*180/pi
 }
}

結果

ユーザ   システム       経過  
  1522.730      1.510   1528.272 

約25分ぐらいですか

次にFortranの場合
コードはこちら

program DistCalc
 implicit none
 real,dimension(1565,5)::Tree
 real,dimension(7146,3)::IncW3
 real,dimension(4952,3)::IncW5
 real,dimension(2108,3)::IncW11
 real pi
 real DistW3(1565*7146),AngW3(1565*7146)
 real DistW5(1565*4952),AngW5(1565*4952)
 real DistW11(1565*2108),AngW11(1565*2108)
 integer i,j 

 open(10, file='Tree')
 open(20, file='IncW3')
 open(30, file='IncW5')
 open(40, file='IncW11')
 read(10,*) ((Tree(i,j),j=1,5),i=1,1565) !iは行、jは列
 read(20,*) ((IncW3(i,j),j=1,3),i=1,7146)
 read(30,*) ((IncW5(i,j),j=1,3),i=1,4952)
 read(40,*) ((IncW11(i,j),j=1,3),i=1,2108)

 do i=1,1565
  do j=1,7146
  DistW3(i*j) = sqrt((Tree(i,4)-IncW3(j,1))**2 + (Tree(i,5)-IncW3(j,2))**2)
  AngW3(i*j) = atan2((Tree(i,5)-IncW3(j,2)),(Tree(i,4)-IncW3(j,1)))
  end do
 end do

 do i=1,1565
  do j=1,4952
  DistW5(i*j) = sqrt((Tree(i,4)-IncW5(j,1))**2 + (Tree(i,5)-IncW5(j,2))**2)
  AngW5(i*j) = atan2((Tree(i,5)-IncW5(j,2)),(Tree(i,4)-IncW5(j,1)))
  end do
 end do

 do i=1,1565
  do j=1,2108
  DistW11(i*j) = sqrt((Tree(i,4)-IncW11(j,1))**2 + (Tree(i,5)-IncW11(j,2))**2)
  AngW11(i*j) = atan2((Tree(i,5)-IncW11(j,2)),(Tree(i,4)-IncW11(j,1)))
  end do
 end do

 pi=4*atan(1.0)

 open(50,file="DistW3")
 do i=1,1565
  do j=1,7146
  write(50,*) Tree(i,4), Tree(i,5), IncW3(j,1), IncW3(j,2),&
   DistW3(i*j), AngW3(i*j), (AngW3(i*j)*180)/pi
  end do
 end do

 open(70,file="DistW5")
 do i=1,1565
  do j=1,4952
  write(70,*) Tree(i,4), Tree(i,5), IncW5(j,1), IncW5(j,2),&
   DistW5(i*j), AngW5(i*j), (AngW5(i*j)*180)/pi
  end do
 end do

 open(80,file="DistW11")
 do i=1,1565
  do j=1,2108
  write(80,*) Tree(i,4), Tree(i,5), IncW11(j,1), IncW11(j,2),&
   DistW11(i*j), AngW11(i*j), (AngW11(i*j)*180)/pi
  end do
 end do

end program DistCalc

で結果はというと、一目瞭然ですね!!

$ f95 DistCalc.f90
$ ./a.out

real 3m20.994s
user 3m14.140s
sys 0m5.330s

今回の場合には、Fortranのコードを書くだけで結構時間がかかった（まだ初心者なので）ので、Rだけで作業しても良かったかもしれません。
ですが、今後のことを考えてもFortranをお勉強したほうが良いですね〜
できればC++もお勉強したいのですが...