|
|
最近在算rotation matrix elemetns相关的，其中有一条是closure的check。具体数学公式之类不详细解释，换成通俗语言问题是这样的:
文件A: 有160个数据，分别是x1(i),y1(i),z1(i)
文件B: 有160个数据，分别是x2(i),y2(i),z2(i)
经过简单运算，二文件数据之差是 x(i),y(i),z(i)
现有S=3(COS(X)^2-1)/2的计算公式，其中COS(X)^2是由x(i),y(i),z(i)运算得来。^是平方符号
现在已经写好的程序是这样的:
[Fortran] 纯文本查看 复制代码 program main
implicit none
integer i,j
integer ios
integer m,n
character(160) :: filename1,filename2,filename3,filename4,tmp
integer id
integer n1,n2
real x1,y1,z1
real x2,y2,z2
real x(10000),y(10000),z(10000)
real v1,v2,v3
real a,b,c
real up,down,coss(10000),s1(10000),s2(10000),s3(10000),avg
integer GetFileN,ii,jj
character(len=10)::status
print*, 'Running...'
filename1='c6_'
filename2='c8_'
do m=1,1 ! number of files to read
write(tmp,*)m
open(10,file=trim(filename1)//trim(adjustl(tmp))
1//'.dat',status='old' )
open(20,file=trim(filename2)//trim(adjustl(tmp))
1//'.dat',status='old' )
do i=1,GetFileN(10)
read(10,*,iostat=ios) n1,x1,y1,z1 !read c6
read(20,*,iostat=ios) n2,x2,y2,z2 !read c8
if (ios /=0) then
exit
endif
x(i)=x2-x1 !c8-c6 x
y(i)=y2-y1 !c8-c6 y
z(i)=z2-z1 !c8-c6 z
c print*,x(i),y(i),z(i)
c write(30,'(3f12.5)')x(i),y(i),z(i)
a=0
b=0
c=1 ! labframe
up=(x(i)*a+y(i)*b+z(i)*c)*(x(i)*a+y(i)*b+z(i)*c)
down=(x(i)*x(i)+y(i)*y(i)+z(i)*z(i))*(a*a+b*b+c*c)
coss(i)=up/down ! cos(i)^2
c print*, coss(i)
s1(i)=(3*coss(i)-1)/2 !n=0 d00
s2(i)=-sqrt(1.5)*sqrt(1-coss(i))*sqrt(coss(i)) ! n=1 d10
s3(i)=sqrt(0.375)*(1-coss(i)) !n=2 d20
c print*, s1(i)
c print*, s2(i)
c print*, s3(i)
enddo
enddo
print*, 'Done.'
close(10)
close(20)
end
integer function GetFileN(iFileUnit)
implicit none
logical , parameter :: b = .True.
integer , intent( IN ) :: iFileUnit
character*(1) :: c
GetFileN = 0
rewind( iFileUnit )
do while (b)
read( iFileUnit , * ,end =999 ,Err = 999 )c
GetFileN = GetFileN + 1
end Do
999 rewind( iFileUnit )
return
end function GetFileN
请忽略n1和n2，不需要他们参与计算。并且对应A和B两种文件其实每种有几万个(所以用了
trim(filename1)//trim(adjustl(tmp))
1//'.dat')
的写法。其他s2和s3是指另外两个需要计算的公式。
因为最后需要对这几万个求出的s1，s2，s3值求平均，用现在的写法计算起来超!慢(这里特别说明一下，因为需要统计上有效，所以实际上还要取不同的间隔。比如文件A1,A2这种间隔为1的，下面还要计算A1,A3这样间隔为2的......直到达到文件总数的一半为止(比方说总共10个文件，实际要重复计算五次，间隔分别是1，2，3，4，5。我之前用过的一种‘土办法’是把每一个n对应的xyz提取出来单独存成另一个文件，用另一个程序计算: 见下面)。所以现在需要超过一天以上的时间完成全部运算!)。我想请教各位有没有别的办法可以实现更快的计算，谢谢!
---土办法---
[Fortran] 纯文本查看 复制代码 print*, 'Please enter the number of files:'
read*,nn
do mm=1,nn
write(tmp,*)mm
open(10,file=trim('S (')//trim(adjustl(tmp))//').dat'
1,status='old')
c open(10,file='nozero3534.dat',status='old')
num=GetFileN(10)
allocate(a(num,3))
allocate(b(num,3))
allocate(c(num))
allocate(x1(num))
allocate(y1(num))
allocate(z1(num))
allocate(x2(num))
allocate(y2(num))
allocate(z2(num))
allocate(up1(num))
allocate(below1(num))
allocate(cos21(num))
allocate(up2(num))
allocate(below2(num))
allocate(cos22(num))
allocate(s(num))
allocate(summ(num))
allocate(avg(num))
do n=1,num
read(10,*,iostat=ios) (a(n,m), m=1,3)
if(ios/=0) then
exit
endif
enddo
do step=1,NUM/2
summ=0
do j=1,NUM-step
B(j,:)=A(step+j,:)
x1(j)=a(j,1)
y1(j)=a(j,2)
z1(j)=a(j,3)
up1(j)=z1(j)*z1(j)
below1(j)=x1(j)*x1(j)+y1(j)*y1(j)+z1(j)*z1(j)
cos21(j)=up1(j)/below1(j)
x2(j)=b(j,1)
y2(j)=b(j,2)
z2(j)=b(j,3)
up2(j)=z2(j)*z2(j)
below2(j)=x2(j)*x2(j)+y2(j)*y2(j)+z2(j)*z2(j)
cos22(j)=up2(j)/below2(j)
s(j)=cos21(j)*cos22(j)
summ=s(j)+summ
enddo
avg=summ(j)/(num-step)
|
|
F 币
贡献
权杖
窥视卡
发表于 2015-3-23 20:19:38
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
楼主
发表于 2015-3-30 08:43:59
