そろそろ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++もお勉強したいのですが...
