一个计算相关的编程思路/方法讨论
最近在算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)运算得来。^是平方符号
现在已经写好的程序是这样的:
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提取出来单独存成另一个文件，用另一个程序计算: 见下面)。所以现在需要超过一天以上的时间完成全部运算!)。我想请教各位有没有别的办法可以实现更快的计算，谢谢!
---土办法---
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)
本帖最后由 li913 于 2015-3-23 22:33 编辑
上面的代码，主题结构没啥可更改的了，适当优化一下。比如
up=(x(i)*a+y(i)*b+z(i)*c)*(x(i)*a+y(i)*b+z(i)*c) 改为：
up=(x(i)*a+y(i)*b+z(i)*c)；up=up*up 。
建议改为动态数组
program main
implicit none
integer i,j
integer ios
integer m,n
character(160) :: filename1,filename2,filename3,filename4,tmp
integer id, n10, n20
integer n1,n2
real x1,y1,z1
real x2,y2,z2
real,allocatable:: x(:),y(:),z(:),s1(:),s2(:),s3(:),coss(:)
real v1,v2,v3
real a,b,c,d
real up,down,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))//'.dat',status='old' )
open(20,file=trim(filename2)//trim(adjustl(tmp))//'.dat',status='old' )
n10=GetFileN(10); n20=GetFileN(20); n10=min(n10,n20)
! 此处分配动态数组
allocate(x(n10),y(n10),z(n10),s1(n10),s2(n10),s3(n10))
a=0; b=0; c=1 ! labframe
d = a*a+b*b+c*c
do i=1,n10
read(10,*) n1,x1,y1,z1 !read c6
read(20,*) n2,x2,y2,z2 !read c8
x(i)=x2-x1 !c8-c6 x
y(i)=y2-y1 !c8-c6 y
z(i)=z2-z1 !c8-c6 z
!up=(x(i)*a+y(i)*b+z(i)*c)*(x(i)*a+y(i)*b+z(i)*c)
up=(x(i)*a+y(i)*b+z(i)*c); up=up*up !优化
! down=(x(i)*x(i)+y(i)*y(i)+z(i)*z(i))*(a*a+b*b+c*c)
down=(x(i)*x(i)+y(i)*y(i)+z(i)*z(i))*d !优化
coss(i)=up/down ! cos(i)^2
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
s2(i)=-sqrt(1.5*(1-coss(i))*coss(i)) !优化
s3(i)=sqrt(0.375)*(1-coss(i)) !n=2 d20
enddo
enddo
print*, 'Done.'
close(10)
close(20)
end li913 发表于 2015-3-23 22:21
上面的代码，主题结构没啥可更改的了，适当优化一下。比如
up=(x(i)*a+y(i)*b+z(i)*c)*(x(i)*a+y(i)*b+z(i) ...
非常感谢!:-handshake
我试了一下你的代码，遇到一个小问题:
http://tinypic.com/view.php?pic=4sno6w&s=8#.VRAqNOEnh_A
对应coss(i)=up/down语句。之前的部分我有逐一检查过。
------
另外就是关于下面的统计部分，您有其他好办法可以让这两个程序写在一起吗? 比如说，对应这样的计算(就是您给优化过的那个)，能不能直接选择file的间隔? 现在是按顺序一个接一个读，A1之后是A2,A3......能不能跳着读? 就是先一个接一个，然后每隔两个，每隔三个，重复同样的计算，一直到全部文件数的一半。我只会下面那个土办法，就是先提取出来，然后用j=1, num-step这样的。这个办法可以用，但是太慢了。我不知道怎样能把它们写在一起。
恳请您的指导!
li913 发表于 2015-3-23 22:21
上面的代码，主题结构没啥可更改的了，适当优化一下。比如
up=(x(i)*a+y(i)*b+z(i)*c)*(x(i)*a+y(i)*b+z(i) ...
为啥不是数据完全读入后，向量化操作呢？边读取边计算的话，编译器没法很好滴优化
lz一点profile的经验也没有？程序热点都找不到的话，等于是无头苍蝇 pasuka 发表于 2015-3-24 08:37
为啥不是数据完全读入后，向量化操作呢？边读取边计算的话，编译器没法很好滴优化
lz一点profile的经验也 ...
您要有好办法就请您用代码介绍介绍，这样我和其他看帖子的人都能尝试运行一下比较哪种更好。我是职业研究数学的，不是职业程序员，学这个就是自己的兴趣。另外麻烦您说话客气点儿。 kif117 发表于 2015-3-24 17:40
您要有好办法就请您用代码介绍介绍，这样我和其他看帖子的人都能尝试运行一下比较哪种更好。我是职业研究 ...
没有人是职业玩 fortran 的。你，我，pasuka 都不是。
楼上的说话一贯就是这样，习惯了就好了。倒是不必太在意了~~
不管是不是职业的，代码里必要的优化还是要做的，是吧？多听听他人的意见，总是能是自己开阔思维。 kif117 发表于 2015-3-24 17:40
您要有好办法就请您用代码介绍介绍，这样我和其他看帖子的人都能尝试运行一下比较哪种更好。我是职业研究 ...
最初搞计算机的可都是学数学出身的，国外的图灵、冯诺依曼不提，国内的103机还是数学家华罗庚在中科院数学所组织的夏培肃老师等人一起开发的
lz这番话把学数学和写代码的这么奚落一番，也真够任性的 fcode 发表于 2015-3-24 21:51
没有人是职业玩 fortran 的。你，我，pasuka 都不是。
楼上的说话一贯就是这样，习惯了就好了。倒是不必 ...
我就是不喜欢上来就摆出教育人的口气说‘你该多看看书’，‘你应该怎样怎样’，‘你怎么样怎么样’，‘你没经验’之类不负责任的话。别人提一个问题，如果你知道又愿意，就请告诉答案。不愿意的话不回答就是了。假设一个学生问我怎么解一道题，他肯定不愿意听见‘你回去看看课本就知道了’这种话，他肯定希望我一步一步解出来给他看解法。公共网路，也应该抱持礼貌客气，您说是吧。 pasuka 发表于 2015-3-24 21:58
最初搞计算机的可都是学数学出身的，国外的图灵、冯诺依曼不提，国内的103机还是数学家华罗庚在中科院数 ...
中国国内我不了解，但我生活的国家不碰编程的纯数学家可是大有人在。有人计算机都不碰，手写论文，学生打进电脑。学习过程中，导师一概要求手脑计算，不得用电脑。我中文不好，但是仔细看一看发言，没觉得哪儿‘任性’。我表达的是: 我是一个专门学习数学知识的人，不擅长编程，也不太花时间去深入钻研编程方法。当然，如果谁有高明的code，还是乐意看一看学一学，可是请您别擅自评论我。倒是您，上来就说我们是‘无头苍蝇’，真的，不知道是谁奚落谁。要是像fcode同学说的那样，您平时就这么说话，我想我们还是别交流的好些。我习惯每天请啊您啊谢谢啊这样的讲话，习惯问题。 kif117 发表于 2015-3-24 22:32
我就是不喜欢上来就摆出教育人的口气说‘你该多看看书’，‘你应该怎样怎样’，‘你怎么样怎么样’，‘你 ...
是的，学生们肯定不愿意听到“回去看看课本”这种回答。
但也就那种就在教科书上明确写着的东西，却有人来问。或者明明是需要整个章节来讲解的知识点，非要让我来给他解释，那我也只能呵呵了......（不指你）
我曾接触过一些人，可能他们就是不喜欢太客气的说话。比如有一个“老师”，就是不允许我叫他老师，并说：“只有你的博硕士导师可以这样称呼”。或者就是不让我称呼他为“您”，而非要我说“你”。
还有一个老外，也很有意思，动辄骂我“stupid”，在我看来很不礼貌。但是他却说：“如果我也犯了这样低级的错误，我很乐意被指责为 stupid。”
尽管如此，上面提到的两位，我还是很尊重他们，并且也发现他们挺可爱的。
最后，确实在编程方面，大家都是业余，这个网站上貌似没有纯学计算机编程的。不过，既然开始写代码了，就把他尽量写完美。
一段优美的代码，不止是可以更准确更快速地得到你的结果，更可以提高你将来再次使用它的频率和效率。对将来你或别人再次阅读它，也会有帮助。（其实说白了，都是方便你自己）