Fortran Coder
标题: 一个计算相关的编程思路/方法讨论 [打印本页]
作者: kif117 时间: 2015-3-23 20:19
标题: 一个计算相关的编程思路/方法讨论
最近在算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)
作者: li913 时间: 2015-3-23 22:21
本帖最后由 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 。
建议改为动态数组
[Fortran] 纯文本查看 复制代码
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
作者: kif117 时间: 2015-3-23 23:28
非常感谢!
我试了一下你的代码，遇到一个小问题:
对应coss(i)=up/down语句。之前的部分我有逐一检查过。
------
另外就是关于下面的统计部分，您有其他好办法可以让这两个程序写在一起吗? 比如说，对应这样的计算(就是您给优化过的那个)，能不能直接选择file的间隔? 现在是按顺序一个接一个读，A1之后是A2,A3......能不能跳着读? 就是先一个接一个，然后每隔两个，每隔三个，重复同样的计算，一直到全部文件数的一半。我只会下面那个土办法，就是先提取出来，然后用j=1, num-step这样的。这个办法可以用，但是太慢了。我不知道怎样能把它们写在一起。
恳请您的指导!
作者: pasuka 时间: 2015-3-24 08:37
为啥不是数据完全读入后，向量化操作呢？边读取边计算的话，编译器没法很好滴优化
lz一点profile的经验也没有？程序热点都找不到的话，等于是无头苍蝇
作者: kif117 时间: 2015-3-24 17:40
您要有好办法就请您用代码介绍介绍，这样我和其他看帖子的人都能尝试运行一下比较哪种更好。我是职业研究数学的，不是职业程序员，学这个就是自己的兴趣。另外麻烦您说话客气点儿。
作者: fcode 时间: 2015-3-24 21:51
没有人是职业玩 fortran 的。你，我，pasuka 都不是。
楼上的说话一贯就是这样，习惯了就好了。倒是不必太在意了~~
不管是不是职业的，代码里必要的优化还是要做的，是吧？多听听他人的意见，总是能是自己开阔思维。
作者: pasuka 时间: 2015-3-24 21:58
最初搞计算机的可都是学数学出身的，国外的图灵、冯诺依曼不提，国内的103机还是数学家华罗庚在中科院数学所组织的夏培肃老师等人一起开发的
lz这番话把学数学和写代码的这么奚落一番，也真够任性的
作者: kif117 时间: 2015-3-24 22:32
我就是不喜欢上来就摆出教育人的口气说‘你该多看看书’，‘你应该怎样怎样’，‘你怎么样怎么样’，‘你没经验’之类不负责任的话。别人提一个问题，如果你知道又愿意，就请告诉答案。不愿意的话不回答就是了。假设一个学生问我怎么解一道题，他肯定不愿意听见‘你回去看看课本就知道了’这种话，他肯定希望我一步一步解出来给他看解法。公共网路，也应该抱持礼貌客气，您说是吧。
作者: kif117 时间: 2015-3-24 22:41
中国国内我不了解，但我生活的国家不碰编程的纯数学家可是大有人在。有人计算机都不碰，手写论文，学生打进电脑。学习过程中，导师一概要求手脑计算，不得用电脑。我中文不好，但是仔细看一看发言，没觉得哪儿‘任性’。我表达的是: 我是一个专门学习数学知识的人，不擅长编程，也不太花时间去深入钻研编程方法。当然，如果谁有高明的code，还是乐意看一看学一学，可是请您别擅自评论我。倒是您，上来就说我们是‘无头苍蝇’，真的，不知道是谁奚落谁。要是像fcode同学说的那样，您平时就这么说话，我想我们还是别交流的好些。我习惯每天请啊您啊谢谢啊这样的讲话，习惯问题。
作者: fcode 时间: 2015-3-25 06:56
是的，学生们肯定不愿意听到“回去看看课本”这种回答。
但也就那种就在教科书上明确写着的东西，却有人来问。或者明明是需要整个章节来讲解的知识点，非要让我来给他解释，那我也只能呵呵了......（不指你）
我曾接触过一些人，可能他们就是不喜欢太客气的说话。比如有一个“老师”，就是不允许我叫他老师，并说：“只有你的博硕士导师可以这样称呼”。或者就是不让我称呼他为“您”，而非要我说“你”。
还有一个老外，也很有意思，动辄骂我“stupid”，在我看来很不礼貌。但是他却说：“如果我也犯了这样低级的错误，我很乐意被指责为 stupid。”
尽管如此，上面提到的两位，我还是很尊重他们，并且也发现他们挺可爱的。
最后，确实在编程方面，大家都是业余，这个网站上貌似没有纯学计算机编程的。不过，既然开始写代码了，就把他尽量写完美。
一段优美的代码，不止是可以更准确更快速地得到你的结果，更可以提高你将来再次使用它的频率和效率。对将来你或别人再次阅读它，也会有帮助。（其实说白了，都是方便你自己）
作者: kif117 时间: 2015-3-25 19:22
不能把缺乏教养和礼貌当成个性，我觉得您认识的那个‘老外’是非常失礼的，在我的国家，这种人会被学生投诉‘歧视’而丢掉工作。他自己喜欢，不代表别人喜欢，也不能作为常识性的标准。我欣赏的中国作家的一位世界级作家，他有一次谈论网络，他说话像写小说一样犀利，他大意说‘人一上网，就脸都不要了’。早些年电脑没有这样普及的时候，bbs不是这样的，绝大多数人都能使用礼貌用语交谈，并且讨论多半都有趣和愉快。现在网路上的中文越来越糟糕了，很多情况词的用法都很古怪，完全读不懂究竟想说什么。我希望在这个论坛上的各位都能尽量使用优美的中文，正确地表达自己的想法。
不知道您有没有教学经验(我自己做过TA)，就算教科书上有，不等于人人可以读。比如一个化学系的人来问问题，我肯定不会建议去读数学书。更不会像某些人那样，‘你一点数学知识都没有吧? Gamma Funtion都不知道还念什么博士? 你是笨蛋菜鸟吗!’就算是再简单的问题，我都只能这样说，‘这个知识点您以前可能不太注意，现在给您讲解一下。’如果我不愿意讲，我会直接拒绝。
---
有没有职业编程的人呢? 我回去仔细想了一下，发现还真有。我的一个同学是制作游戏的，她每天的工作就是编写调试程序。当然，这个‘职业’，是指靠此谋生的意思。如果让她自己说，很可能回答职业是‘饲养宠物猫’。
作者: fcode 时间: 2015-3-25 19:57
中国自古就是礼仪之邦，其实几千年来，大多数人还是比较讲究礼仪的。日本的文化里也是比较注重教养和礼仪的。
不过，我会乐于接纳其他人的性格特点。
我想，当面骂我的人，总比当面对我笑呵呵而背地里骂我的好。
我的教学经验不多，做过一段时间的售后服务支持。就这种工作性质来说，我肯定不会真正对“上帝”发火。但坦白的说，面对一些基础比较薄弱的客户，虽然嘴上不说，心里还是想骂几句的。
---
当然，有职业的程序员，而且很多。但是，我说的是玩 fortran 的，这么多年了，还真没遇到过。所以，这个论坛上，大家都是编程的“外行”。
作者: andy8496 时间: 2015-3-25 20:03
我倒是很喜欢pasuka这种风格
就是不知道ta是男是女？
每次被虐过后自己再恶补一番，总有技能唰唰唰往上飙的感觉。
最怕的是一些高手见了简单问题不理不睬的飘过。
感谢各位高手手把手的教导！感谢fcode伴我一路成长！
这是我见过风气最好的论坛，没有之一。
作者: kif117 时间: 2015-3-25 21:17
我想这是文化差异问题，在我生活的国家，就算学习好别人也不会觉得你有权力对学习不好的人做负面评价。对博士生来说，(我们组的要求)是先做好人，再做学问。无论跟谁都用sie(您)讲话，这样即使小朋友听了，也会自然学习优美的语言，将它变成一种习惯。我们的社会通常比较单纯，也就是中国人说的‘傻老外’，做学问的人通常都很和气。就算学生问了一个真正的蠢问题，也不会觉得恼火，你可以反过来想，现在他总算知道正确答案了。要知道别人可能在其他方面很优秀，有游泳的世界冠军，短跑的世界冠军，但这些人在您和我的专业范围内几乎是0基础。倘若他们愿意加入，我认为我们应该宽容地对待。
现在像您这样能够认真思考并且在(中文为主的)论坛上发表意见的人不太多了，要知道，我第一次收到’不明觉厉‘的回答时候甚至跑去查字典! 后来发现这在中文论坛是’基本常识‘。
方便的话，我能麻烦您看一下这个程序吗 (就是主题里提到的’土办法‘)? 对于9300个数据进行计算，花费1089s。我希望能更快的运算，但没有(能够)使得它变好的任何经验。目前来说我可以完成一项普通程度的编程工作，但不能用最短时间完成，而这也是我很感兴趣的部分。
[Fortran] 纯文本查看 复制代码
program main
implicit none
integer,parameter::file_number= 1
character(100)::tmp
character(160)::fileN
character(100)::filename(file_number)
integer::i,j,k,m,n,NUM,step,ios,ntotal
integer::p,q,GetFileN
real,pointer::A(:,:),B(:,:),c(:)
real,pointer::x1(:)
real,pointer::y1(:)
real,pointer::z1(:)
real,pointer::x2(:)
real,pointer::y2(:)
real,pointer::z2(:)
real,pointer::up(:)
real,pointer::below(:)
real,pointer::cos2(:)
real,pointer::s(:)
real,pointer::summ(:)
real,pointer::avg(:)
real*8 time
integer*4 time0, time1, dtime
c Start time
call system_clock(time0)
open(10,file='testdata.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(up(num))
allocate(below(num))
allocate(cos2(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
c write(6,'(f8.3,f8.3,f8.3)')(a(n,m), m=1,3)
enddo
Print*, 'Running...'
do step=1,NUM/2
summ=0
do j=1,NUM-step
c write(tmp,"(i4.4)")step
B(j,:)=A(step+j,:)
c print*, A(j,:)
c print*, B(j,:)
x1(j)=a(j,1)
y1(j)=a(j,2)
z1(j)=a(j,3)
x2(j)=b(j,1)
y2(j)=b(j,2)
z2(j)=b(j,3)
up(j)=(x1(j)*x2(j)+y1(j)*y2(j)+z1(j)*z2(j))**2
below(j)=(x1(j)*x1(j)+y1(j)*y1(j)+z1(j)*z1(j))*(x2(j)*x2(j)
1+y2(j)*y2(j)+z2(j)*z2(j))
cos2(j)=up(j)/below(j)
s(j)=(3*cos2(j)-1)/2
summ=s(j)+summ
enddo
c print*, summ(j)
avg=summ(j)/(num-step)
c print*, (num-step)
c print*, avg(j)
open(20,file='testresult.dat',status='unknown')
write(20,222)avg(j)
enddo ! enddo for 'tmp'=>step
111 format (3f8.3)
222 format (f12.6)
close(10)
deallocate(A,B,C)
deallocate(x1,y1,z1,x2,y2,z2,up,below,cos2,s,summ,avg)
c Finish time
call system_clock(time1, dtime)
c Output time
time = 1d0*(time1-time0)/dtime
write(*,"(a7,f16.7)")"Time = ",time
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
数据是这样的:
0 0 1
21.086 0.082 21.079
20.531 -0.909 20.161
19.106 -0.566 19.855
18.938 -0.763 21.644
19.99 -2.14 22.666
...
(total=9300)
感谢您花费时间，并且感谢您的帮忙!
作者: fcode 时间: 2015-3-25 21:23
我建议你把 9300 个数据传上来。
（编辑框里如果有字数限制，可以用附件上传）
作者: kif117 时间: 2015-3-25 21:33
以下是全部的9300个数据，再次感谢您的帮助!
---
第一次学会发送‘附件’!
-
-
testdata.TXT
189.07 KB, 下载次数: 9
作者: kif117 时间: 2015-3-25 21:35
全部数据已经用附件上传，经过我个人测试是可以正确下载的。
再次感谢您的帮忙。
作者: fcode 时间: 2015-3-25 22:06
Running...
Time = 7.8090000
在我的计算机上，修改后的代码计算 9300 个数用了 7 秒，不知道你是否可以接受？
如果还要继续优化，恐怕就要从算法上下功夫了。看起来有点像“自相关”？
[Fortran] 纯文本查看 复制代码
Program main
Implicit None
Integer, Parameter :: file_number = 1
Character (100) :: tmp
Character (160) :: filen
Character (100) :: filename(file_number)
Integer :: i, j, k, m, n, num, step, ios, ntotal
Integer :: p, q, GetFileN
Real, Pointer :: a(:, :), b(:, :), c(:)
Real, Pointer :: x1(:)
Real, Pointer :: y1(:)
Real, Pointer :: z1(:)
Real, Pointer :: x2(:)
Real, Pointer :: y2(:)
Real, Pointer :: z2(:)
!Real, Pointer :: up(:)
!Real, Pointer :: below(:)
!Real, Pointer :: cos2(:)
!Real, Pointer :: s(:)
!Real, Pointer :: summ(:)
Real :: summ , avg , s , up , below , cos2
!Real, Pointer :: avg(:)
Real *8 time
Integer *4 time0, time1, dtime
Call system_clock(time0)
Open (10, File='testdata.dat', Status='old')
Open (20, File='testresult.dat', Status='unknown')
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 (up(num))
!Allocate (below(num))
!Allocate (cos2(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) Exit
End Do
Print *, 'Running...'
x1(:) = a(:,1)
y1(:) = a(:,2)
z1(:) = a(:,3)
Do step = 1, num/2
summ = 0.0
avg = 0.0
!write(*,*) step
Do j = 1, num - step
b(j, :) = a(step+j, :)
x2(j) = b(j, 1)
y2(j) = b(j, 2)
z2(j) = b(j, 3)
up = (x1(j)*x2(j)+y1(j)*y2(j)+z1(j)*z2(j))**2
below = (x1(j)*x1(j)+y1(j)*y1(j)+z1(j)*z1(j))*(x2(j)*x2(j)+y2(j)*y2(j)+z2(j)*z2(j))
cos2 = up/below
s = (3*cos2-1)/2
summ = s + summ
End Do
avg = summ/(num-step)
Write (20, 222) avg
End Do
Close (10)
Close (12)
Deallocate (a, b, c)
Deallocate (x1, y1, z1, x2, y2, z2 ) !, up, below, cos2 , s , summ, avg)
Call system_clock(time1, dtime)
time = 1D0*(time1-time0)/dtime
Write (*, '(a7,f16.7)') 'Time = ', time
111 Format (3F8.3)
222 Format (F12.6)
End Program main
Integer function GetFileN(iFileUnit)
implicit none
logical , parameter :: b = .True.
integer , intent( IN ) :: iFileUnit
character(len=1) :: c
GetFileN = 0
Rewind( iFileUnit )
Do while (b)
Read( iFileUnit , * ,end =999 ,Err = 999 )c
GetFileN = GetFileN + 1
End Do
999 Rewind( iFileUnit )
End function GetFileN
作者: kif117 时间: 2015-3-25 22:23
试验以后难以表达我的喜悦之感!!! 现在我可以用fs为单位进行计算了(对精确的追求方面有些神经质)。您说得对，这个是correlation function (绝大多数专业词我不会用中文表达，非常抱歉)相关的计算。就是Satchler那本Anguler Momentum147页提到的closure部分的验证。这个角度cos(theta)是stochastic quantity里面的一个时间相关的参数(我表达的不好，可能词不达意)。为了统计上的精确，不但时间间隔取1，还要取2，取3.....假设时间为10，那么要取到至少间隔5; 就是用第一个数据和第六个数据计算，第二个数据和第七个数据计算.....以此类推。就结果本身来说其实没有特别大的差别，但如果plot出来就会发现如果不这样图像会出现锯齿，而不是smooth的。
---
我能不能请教您，我其实看不太出来究竟是哪里的修改起了关健性的作用。为什么计算时间会相差上百倍? 能不能请您简单地谈一谈，根据经验，如何能在做这类步骤很多的重复的计算时优化速度。非常感谢!!!
作者: fcode 时间: 2015-3-25 22:32
我想你之前的代码可能最严重的错误就是
summ=s(j)+summ
这一句了。
这一句本身就在两重循环里，被执行的次数非常多。而 summ 又是一个数组，相当于三重循环了。
而实际上，summ 只需要一个数就行了。
当然，你把这一句改成
summ(j)=s(j)+summ(j)
也可以。
不过，如果能用单变量，就不要用数组，毕竟内存占用差了 9000 倍。
另外就是
x1(:) = a(:,1)
y1(:) = a(:,2)
z1(:) = a(:,3)
这三句被执行了大概 num/2 * num 次，实际上也只需要执行一次，放到循环最前面就行了。
correlation function 应该就是“相关”了，autocorrelation 在中文里叫自相关。
如果是一维的情况，数据量大一些我会选择在频率域进行。
不过是否适合你的三维情况，还需要推导一下才知道。
作者: kif117 时间: 2015-3-25 23:23
我这里还有一个相关联的很老的程序(可能是我出生以前出版的)。因为不明白什么是‘NUMBER OF STEPS ON THE TAPE’，所以也没有试验过。也许像您这样Fortran经验丰富的人可以看懂。
********************************************************************************
** FICHE F.27. PROGRAM TO COMPUTE TIME CORRELATION FUNCTIONS **
** This FORTRAN code is intended to illustrate points made in the text. **
** To our knowledge it works correctly. However it is the responsibility of **
** the user to test it, if it is to be used in a research application. **
********************************************************************************
********************************************************************************
** FICHE F.27. PROGRAM TO COMPUTE TIME CORRELATION FUNCTIONS **
** This FORTRAN code is intended to illustrate points made in the text. **
** To our knowledge it works correctly. However it is the responsibility of **
** the user to test it, if it is to be used in a research application. **
********************************************************************************
PROGRAM TCORR
COMMON / BLOCK1 / STORX, STORY, STORZ
COMMON / BLOCK2 / VX, VY, VZ
COMMON / BLOCK3 / VACF, ANORM
C *******************************************************************
C ** CALCULATION OF TIME CORRELATION FUNCTIONS. **
C ** **
C ** THIS PROGRAM ANALYZES DATA TO CALCULATE A TIME CORRELATION **
C ** FUNCTION IN ONE SWEEP (LOW MEMORY REQUIREMENT). IN THIS **
C ** EXAMPLE THE VELOCITY AUTO-CORRELATION FUNCTION IS CALCULATED. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** INTEGER N NUMBER OF ATOMS **
C ** INTEGER NSTEP NUMBER OF STEPS ON THE TAPE **
C ** INTEGER IOR INTERVAL FOR TIME ORIGINS **
C ** INTEGER NT CORRELATION LENGTH, INCLUDING T=0 **
C ** INTEGER NTIMOR NUMBER OF TIME ORIGINS **
C ** INTEGER NLABEL LABEL FOR STEP (1,2,3...NSTEP) **
C ** REAL VX(N),VY(N),VZ(N) VELOCITIES **
C ** REAL VACF(NT) THE CORRELATION FUNCTION **
C ** NSTEP AND NT SHOULD BE MULTIPLES OF IOR. **
C ** **
C ** ROUTINES REFERENCED: **
C ** **
C ** SUBROUTINE STORE ( J1 ) **
C ** ROUTINE TO STORE THE DATA FOR CORRELATION **
C ** SUBROUTINE CORR ( J1, J2, IT ) **
C ** ROUTINE TO CORRELATE THE STORED TIME ORIGINS **
C ** **
C ** USAGE: **
C ** **
C ** DATA IN FILE DFILE ON FORTRAN UNIT DUNIT. **
C ** RESULTS IN FILE RFILE ON FORTRAN UNIT RUNIT. **
C *******************************************************************
INTEGER N, NSTEP, IOR, NT, NDIM, DUNIT, RUNIT, NTIMOR
INTEGER FULLUP
PARAMETER ( N = 9300, NSTEP = 9299, IOR = 1, NT = 9300 )
PARAMETER ( DUNIT = 10, RUNIT = 11 )
PARAMETER ( NDIM = NT / IOR + 1, NTIMOR = NSTEP / IOR )
PARAMETER ( FULLUP = NDIM - 1 )
REAL VX(N), VY(N), VZ(N)
REAL STORX(NDIM,N), STORY(NDIM,N), STORZ(NDIM,N)
REAL VACF(NT), ANORM(NT)
INTEGER S(NTIMOR), TM(NTIMOR)
INTEGER TS, TSS, L, NINCOR, K, JA, IB, IN, IA, JO, I
INTEGER NLABEL
CHARACTER DFILE * 30
CHARACTER RFILE * 30
C *******************************************************************
WRITE(*,'('' **** PROGRAM TCORR **** '')')
WRITE(*,'('' CALCULATION OF TIME CORRELATION FUNCTIONS '')')
C ** READ IN FILE NAMES **
WRITE(*,'('' ENTER DATA FILE NAME '')')
READ (*,'(A)') DFILE
WRITE (*,'('' ENTER RESULTS FILE NAME '')')
READ (*,'(A)') RFILE
C ** INITIALIZE COUNTERS **
NINCOR = FULLUP
JA = 1
IA = 1
IB = 1
C ** ZERO ARRAYS **
DO 5 I = 1, NT
VACF(I) = 0.0
ANORM(I) = 0.0
5 CONTINUE
C ** OPEN DATA FILE AND RESULTS FILE **
OPEN ( UNIT = DUNIT, FILE = DFILE, ACCESS = 'SEQUENTIAL',
: STATUS = 'OLD', FORM = 'UNFORMATTED' )
OPEN ( UNIT = RUNIT, FILE = RFILE, STATUS = 'UNKNOWN' )
C ** CALCULATION BEGINS **
DO 40 L = 1, NTIMOR
JA = JA + 1
S(L) = JA - 1
READ ( DUNIT ) NLABEL, VX, VY, VZ
TM(L) = NLABEL
C ** STORE STEP AS A TIME ORIGIN **
CALL STORE ( JA )
C ** CORRELATE THE ORIGINS IN STORE **
DO 10 IN = IA, L
TSS = TM(L) - TM(IN)
TS = TSS + 1
JO = S(IN) + 1
CALL CORR ( JO, JA, TS )
10 CONTINUE
C ** READ IN DATA BETWEEN TIME ORIGINS. THIS CAN **
C ** BE CONVENIENTLY STORED IN ELEMENT 1 OF THE **
C ** ARRAYS STORX ETC. AND CAN THEN BE CORRELATED **
C ** WITH THE TIME ORIGINS. **
DO 30 K = 1, IOR - 1
READ ( DUNIT ) NLABEL, VX, VY, VZ
CALL STORE ( 1 )
DO 20 IN = IA, L
TSS = NLABEL - TM(IN)
TS = TSS + 1
JO = S(IN) + 1
CALL CORR ( JO, 1, TS )
20 CONTINUE
30 CONTINUE
IF ( L .GE. FULLUP ) THEN
IF ( L .EQ. NINCOR ) THEN
NINCOR = NINCOR + FULLUP
JA = 1
ENDIF
IA = IA + 1
ENDIF
40 CONTINUE
CLOSE ( UNIT = DUNIT )
C ** NORMALISE CORRELATION FUNCTIONS **
VACF(1) = VACF(1) / ANORM(1) / REAL ( N )
DO 50 I = 2, NT
VACF(I) = VACF(I) / ANORM(I) / REAL ( N ) / VACF(1)
50 CONTINUE
WRITE ( RUNIT, '('' VELOCITY ACF '')')
WRITE ( RUNIT, '(I6,E15.6)') ( I, VACF(I), I = 1, NT )
CLOSE ( RUNIT )
STOP
END
SUBROUTINE STORE ( J1 )
COMMON/ BLOCK1 / STORX, STORY, STORZ
COMMON/ BLOCK2 / VX, VY, VZ
C *******************************************************************
C ** SUBROUTINE TO STORE TIME ORIGINS **
C *******************************************************************
INTEGER J1
INTEGER N, NT, IOR, NDIM
PARAMETER ( N = 256, NT = 200, IOR = 4 )
PARAMETER ( NDIM = NT / IOR + 1 )
REAL STORX(NDIM,N), STORY(NDIM,N), STORZ(NDIM,N)
REAL VX(N), VY(N), VZ(N)
INTEGER I
DO 10 I = 1, N
STORX(J1,I) = VX(I)
STORY(J1,I) = VY(I)
STORZ(J1,I) = VZ(I)
10 CONTINUE
RETURN
END
SUBROUTINE CORR ( J1, J2, IT )
COMMON/ BLOCK1 / STORX, STORY, STORZ
COMMON/ BLOCK3 / VACF, ANORM
C *******************************************************************
C ** SUBROUTINE TO CORRELATE TIME ORIGINS **
C *******************************************************************
INTEGER J1, J2, IT
INTEGER N, NT, IOR, NDIM
PARAMETER ( N = 256, NT = 200, IOR = 4 )
PARAMETER ( NDIM = NT / IOR + 1 )
REAL STORX(NDIM,N), STORY(NDIM,N), STORZ(NDIM,N)
REAL VACF(NT), ANORM(NT)
INTEGER I
C *******************************************************************
DO 10 I = 1, N
VACF(IT) = VACF(IT) + STORX(J1,I) * STORX(J2,I)
: + STORY(J1,I) * STORY(J2,I)
: + STORZ(J1,I) * STORZ(J2,I)
10 CONTINUE
ANORM(IT) = ANORM(IT) + 1.0
RETURN
END
作者: kif117 时间: 2015-3-26 03:35
另外，我参考了您数据分类问题求助的帖子，改写了一下提取坐标的程序(之前我是用IF判断ID然后取它的对应坐标，是比较慢的方法。)，但是遇到一点问题:
[Fortran] 纯文本查看 复制代码
Program main
Implicit None
Character :: Filename
Character :: Filename1
Integer :: I, J, num
real :: A(4), B(4), C(4)
Open (10, File='IDlist.dat')
num=160
Do J = 1, num ! number of files
Read (10, *) Filename1
Open (J, File='c6c8_ID_'//Filename1//'.dat') ! save result here
write (J,'(3f8.3)') 0.0,0.0,1.0
End Do
Open (20, File='filename.dat')
Do I = 1, num
Read (20, *) Filename
Open (30, File='c6_'//Filename//'.dat')
Open (40, File='c8_'//Filename//'.dat')
Do J = 1, 2
Read (30, *) A(:)
Read (40, *) B(:)
C(:)=B(:)-A(:)
Write (A(1), '(3f8.3)') C(2:4)
End Do
Close (30)
Close (40)
End Do
Close (20)
Do J = 1, num
Close (J)
End Do
Close (10)
End Program
这是我参考后写的程序。对应每个单一文件，总共有160个ID以及它们的坐标。问题是: 对于每个Cn_n各有92000个! 假设它们分别位于D:\02.Fortran1 和D:\02.Fortran2 文件夹下，如何在不同的文件夹下打开相应的文件? 尝试了‘D:\02.Fortran1’//trim(Cn_)//Filename这样的写法，但是不成功。另外，如果我将全部两组各92000个文件放在同一个文件夹下运行，就会出现‘orrt1: The process cannot access the file because it is being used by another pocess:...’这样从未遇到过的错误。
之前您帮我优化的部分，现在写成这样做批量计算:
[Fortran] 纯文本查看 复制代码
Call system_clock(time0)
tn=2 ! here insert the total files you want to run
Do i = 1,tn
write(tmp,*)i
open(10,file=trim('c6c8_ID_')//trim(adjustl(tmp))
1//'.dat',status='old' )
open(20,file=trim('R_')//trim(adjustl(tmp))
1//'.dat',status='unknown' )
num = GetFileN(10)
End Do
Allocate (a(num,3))
Allocate (b(num,3))
Allocate (c(num))
Allocate (x1(num))
我不知道能不能把这两个程序写在一起(当然我希望能够)，我通常只能每个步骤(提取坐标，再计算)分别运行。不知道高阶的编程者是倾向于拥有全部功能的程序，还是为了检查方便而分开写?
谢谢您的帮助! Thank you very much!
作者: kif117 时间: 2015-3-26 23:11
补充: 当我加入批量计算
tn=2 ! //here insert the total files you want to run
Do i = 1,tn
write(tmp,*)i
open(10,file=trim('c6c8_ID_')//trim(adjustl(tmp))//'.dat',status='old' )
open(20,file=trim('R_')//trim(adjustl(tmp))//'.dat',status='unknown' )
num = GetFileN(10)
End Do
部分后，每次都只能保存最后一次计算结果，加上position='append'也无法解决问题。
作者: fcode 时间: 2015-3-27 08:30
1. 关于 22 楼的代码，一个非我的专业的，且没有任何注释的代码，我是无法看懂的。最多在语法层次上看懂，但具体有什么物理意义，就不知道了。
2. 尝试用 "D:\02.Fortran1\" // trim(Cn_) // trim(Filename)
3. The process cannot access the file because it is being used by another pocess:... 的错误应该是无法打开或读写文件，因为另一个进程正在占用它。（另一个进程可能也是自身，比如 Open 以后没有关闭就打开相同的文件）
4. 有必要的话，我会把程序写在一起。而不是选择分开写。您可以学习一下新的 Fortran90 语法，例如 Module 这些，会更方便你对一个大程序进行模块划分。
作者: kif117 时间: 2015-3-27 18:20
首先我必须感谢您的帮助，像这样可以请教具体的问题细节，似乎变得对编程更加感兴趣起来!
Program main
Implicit None
Character :: Filename
Character :: Filename1
character(4)::fn2,fn3
Integer :: I, J,num
real :: A(4), B(4), C(4)
num=9 ! number of files you want to run
Open (10, File='IDlist.dat')
Do J = 1, num
Read (10, *) Filename1
Open (J, File='c6c8_ID_'//Filename1//'.dat',status='unknown')
! save result here
write (J,'(3f8.3)') 0.0,0.0,1.0
End Do
fn2='c6_'
fn3='c8_'
Open (20, File='filename.dat')
Do I = 1, num ! number of files you want to run
Read (20, *) Filename
Open (30, File="D:\C6\"//trim(fn2) //Filename//'.dat')
Open (40, File="D:\C8\"//trim(fn3) //Filename//'.dat')
Do J = 1, num ! number of files you want to run
Read (30, *) A(:)
Read (40, *) B(:)
C(:)=B(:)-A(:)
Write (A(1), '(3f8.3)') C(2:4)
End Do
Close (30)
Close (40)
End Do
Close (20)
Do J = 1, num
Close (J)
End Do
Close (10)
End Program
现在的问题是num只能取到9，变成10，就会出现无法打开文件的错误。我试验了改序列号(比如open 1000...)之类的，但是仍然出现同样的错误。遇到这样的问题我就无能为力了，目前我的水平只能改一些compile时遇到的基本错误。能不能请您就这个程序再次指导一下? 这个程序是我参考了您那篇指导关于提取文件中的年份做文件名，并将年份后面的数值写进同对应年份文件的例子后写的。以前遇到这样的问题都是用if语句判断，运算起来自然没有您的方法快。
另外遇到的一件事是，不能用// trim(Cn_) //这样的写法，我必须写成filename=xxx，然后再用trim(filename)的格式。也不明白原因。
献上诚挚的谢意!
作者: 楚香饭 时间: 2015-3-27 19:58
看注释
Program main
Implicit None
Character(len=8) :: filename !// 字符串要给长度
Character(len=8) :: filename1 !// 如果不给长度，默认是1个长度。只能到9，10就不行了
Integer :: i, j, num
Real :: a(4), b(4), c(4)
num=9 ! number of files you want to run
Open (10, File='IDlist.dat')
Do j = 1, num
Read (10, *) filename1
Open (j, File='c6c8_ID_'//filename1//'.dat', Status='unknown')
Write (j, '(3f8.3)') 0.0, 0.0, 1.0
End Do
Open (20, File='filename.dat')
Do i = 1, num ! number of files you want to run
Read (20, *) filename
Open (30, File='D:\C6\c6_'//trim(adjustl(filename))//'.dat') !// filename 最好 adjustl一下，去除左边的空格
Open (40, File='D:\C8\c8_'//trim(adjustl(filename))//'.dat') !// trim 只能去除右边的空格
Do j = 1, num ! number of files you want to run
Read (30, *) a(:)
Read (40, *) b(:)
c(:) = b(:) - a(:)
Write (a(1), '(3f8.3)') c(2:4)
End Do
Close (30)
Close (40)
End Do
Close (20)
Do j = 1, num
Close (j)
End Do
Close (10)
End Program main
作者: kif117 时间: 2015-3-27 20:28
楚香饭 发表于 2015-3-27 19:58
看注释<pre class="brush:fortran;">Program main Implicit None Character(len=8) ...
Character(160) :: filename, filename1
character(160):: fn1, fn2, fn3
Integer :: I, J,num
real :: A(4), B(4), C(4)
num=10 ! number of files you want to run
fn1='c6c8_'
fn2='c6_'
fn3='c8_'
Open (10, File='IDlist.dat')
Do J = 1, num
Read (10, *) Filename1
Open (J, File="D:\C6c8\"
1//trim(adjustl(('C8_')))//trim(adjustl((filename1)))//'.dat',
1status='unknown')
! save result here
write (J,'(3f8.3)') 0.0,0.0,1.0
End Do
Open (20, File='filename.dat')
Do I = 1, num ! number of files you want to run
Read (20, *) Filename
Open (30, File="D:\C6\"
1//trim(adjustl((fn2)))//trim(adjustl((filename)))//'.dat',
1status='unknown')
Open (40, File="D:\C8\"
1//trim(adjustl((fn3)))//trim(adjustl((filename)))//'.dat',
1status='unknown')
Do J = 1, num ! number of files you want to run
...
我按照您的建议改了一下(索性设置成160， 全部使用trin(adjustl(...))的写法)，现在可以运行到10，一旦设置NUM大于10，就会提示'end-of-file during read unit 10.....c6c8_10.dat‘。能再次请教您吗?
作者: 楚香饭 时间: 2015-3-27 20:34
说明你的 c6c8_10.dat 里面没有 10 行呗。
Do J = 1, num 说明至少要读 num 行
作者: chiangtp 时间: 2015-3-29 01:38
接力14-20樓的討論, 請kif117參考
test.f90
(1.98 KB, 下载次数: 4)
作者: kif117 时间: 2015-3-29 04:06
我单纯地汇报一下在我的计算机上运行的时间:
fcode帮助修改我自己写的‘土程序’的运行时间: 2.87s
您帮我写的例题程序的运行时间: 1.14s
我想请教一下，语句REWIND()，如果我需要修改成批量打开文件TRIM(ADJUSTL(testdata))'.dat'这样，这里要怎么处理? 它定义是'statement positions the file associated with the specified unit'...批量打开的时候，比如我写do i=1,num，是rewind(unit=i)这样吗?
以及请您看一下26楼里的这句‘ Write (a(1), '(3f8.3)') c(2:4)’请问a(1)是指什么? 我没有太看懂这里。所以越改越错.....
谢谢!
------
在遇到同学们以后这几日我每天都花费一两个小时在编程上，感觉像玩LEGO积木那样‘越玩越想玩’。可见，在有指导的情况下学习是很愉快的! 非常感谢大家诚恳的帮助!
作者: fcode 时间: 2015-3-29 07:59
修改为 TRIM(ADJUSTL(testdata)) // '.dat'
rewind(文件通道号) 就可以了，文件通道号需要事先 Open 语句里指定。
Write (a(1), '(3f8.3)') c(2:4) 写法不对，是你先在22楼给出来的。
作者: kif117 时间: 2015-3-30 03:40
Program Www_fcode_cn
Implicit None
Character *5 Filename
Character *4 Filename1
Integer :: I, J
Integer :: A(6)
Open (10, File='D:\a\year.txt')
Do J = 1961, 2013
Read (10, *) Filename1
Open (J, File='D:\c\'//Filename1//'.txt') !// 11 改为 j
End Do
Open (12, File='D:\a\tem.txt')
Do I = 1, 525
Read (12, *) Filename
Open (13, File='D:\a\'//Filename//'.txt')
Do J = 1961, 2013
Read (13, *) A(:)
Write (A(5), *) A(:)
End Do
Close (13)
End Do
Close (12)
Do J = 1961, 2013
Close (J)
End Do
Close (10)
End Program Www_fcode_cn
这个是您帮助另一个同学改写的，我拿来参考(我们遇到的问题是几乎一样的，她是提取特定的年份，我是提取特定的坐标作为result文件的名字):
Program main
Implicit None
Character(160) :: Filename
Character(160) :: fn1,fn2,fn3,tmp
Integer :: I, J,num
real :: A(4), B(4), C(4)
integer ios
fn1='c6c8_ID_'
fn2='c6_'
fn3='c8_'
num=160 ! total number of result files
print*, 'Running'
do J = 1, num
write (tmp,*)J
C print*, J
open (j,File="D:\c6c8\"//trim(adjustl(fn1))//trim(adjustl(tmp))//'.dat',status='unknown') ! save result here
write (j,'(3f8.3)') 0.0,0.0,1.0
enddo
open (20, File='IDlist.dat',status='old') ! 里面存放的是10的倍数
do I = 1, 9200
read (20, *,iostat=ios) Filename
if (ios /=0) then
exit
endif
c print*, filename
open (200,File="D:\C6\"//trim(adjustl(fn2)) //trim(adjustl(filename))//'.dat')
open (300,File="D:\C8\"//trim(adjustl(fn3)) //trim(adjustl(filename))//'.dat')
Do J = 1, num
Read (200, *) A(:) ! first row=ID number， then x,y,z
Read (300, *) B(:) ! first row=ID number
C(:)=B(:)-A(:)
c print*, c(:)
Write (c(1), '(3f8.3)') C(2:4) ! first row=ID number
c End Do
Close (300)
Close (200)
enddo
enddo
Close (20)
Do J = 1, num
Close (J)
End Do
End Program
从周五开始我已经连续改了三天，不知道是什么错误。非常想学习如何不使用(if ID=特定数)的判定语句的情况下达成同样result的程序的写法，特别向您请教程序问题的所在!
作者: fcode 时间: 2015-3-30 08:43
我建议你单独开一个帖子来讨论这个问题。
给出你的输入文件（范例，比如里面只有两行），然后你的输出文件应该是怎么样的？重点说明你的需求。表达明确。
我很努力的从你的代码里看出你的意图，但我感觉比较昏乱，而且可能存在逻辑错误。所以，我对你的需求还不了解。
作者: kif117 时间: 2015-3-30 19:01
好的，我另外发了一个新的帖子，麻烦您指导一下。
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) |
Powered by Discuz! X3.2 |