一个计算相关的编程思路/方法讨论

kif117

不能把缺乏教养和礼貌当成个性，我觉得您认识的那个‘老外’是非常失礼的，在我的国家，这种人会被学生投诉‘歧视’而丢掉工作。他自己喜欢，不代表别人喜欢，也不能作为常识性的标准。我欣赏的中国作家的一位世界级作家，他有一次谈论网络，他说话像写小说一样犀利，他大意说‘人一上网，就脸都不要了’。早些年电脑没有这样普及的时候，bbs不是这样的，绝大多数人都能使用礼貌用语交谈，并且讨论多半都有趣和愉快。现在网路上的中文越来越糟糕了，很多情况词的用法都很古怪，完全读不懂究竟想说什么。我希望在这个论坛上的各位都能尽量使用优美的中文，正确地表达自己的想法。

不知道您有没有教学经验(我自己做过TA)，就算教科书上有，不等于人人可以读。比如一个化学系的人来问问题，我肯定不会建议去读数学书。更不会像某些人那样，‘你一点数学知识都没有吧? Gamma Funtion都不知道还念什么博士? 你是笨蛋菜鸟吗!’就算是再简单的问题，我都只能这样说，‘这个知识点您以前可能不太注意，现在给您讲解一下。’如果我不愿意讲，我会直接拒绝。

---
有没有职业编程的人呢? 我回去仔细想了一下，发现还真有。我的一个同学是制作游戏的，她每天的工作就是编写调试程序。当然，这个‘职业’，是指靠此谋生的意思。如果让她自己说，很可能回答职业是‘饲养宠物猫’。

fcode

中国自古就是礼仪之邦，其实几千年来，大多数人还是比较讲究礼仪的。日本的文化里也是比较注重教养和礼仪的。

不过，我会乐于接纳其他人的性格特点。

我想，当面骂我的人，总比当面对我笑呵呵而背地里骂我的好。

我的教学经验不多，做过一段时间的售后服务支持。就这种工作性质来说，我肯定不会真正对“上帝”发火。但坦白的说，面对一些基础比较薄弱的客户，虽然嘴上不说，心里还是想骂几句的。

---
当然，有职业的程序员，而且很多。但是，我说的是玩 fortran 的，这么多年了，还真没遇到过。所以，这个论坛上，大家都是编程的“外行”。

andy8496

我倒是很喜欢pasuka这种风格
就是不知道ta是男是女？
每次被虐过后自己再恶补一番，总有技能唰唰唰往上飙的感觉。
最怕的是一些高手见了简单问题不理不睬的飘过。
感谢各位高手手把手的教导！感谢fcode伴我一路成长！
这是我见过风气最好的论坛，没有之一。

kif117

fcode 发表于 2015-3-25 19:57
中国自古就是礼仪之邦，其实几千年来，大多数人还是比较讲究礼仪的。日本的文化里也是比较注重教养和礼仪的 ...

我想这是文化差异问题，在我生活的国家，就算学习好别人也不会觉得你有权力对学习不好的人做负面评价。对博士生来说，(我们组的要求)是先做好人，再做学问。无论跟谁都用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

我建议你把 9300 个数据传上来。
（编辑框里如果有字数限制，可以用附件上传）

kif117

以下是全部的9300个数据，再次感谢您的帮助!

---
第一次学会发送‘附件’!

kif117

fcode 发表于 2015-3-25 21:23
我建议你把 9300 个数据传上来。
（编辑框里如果有字数限制，可以用附件上传） ...

全部数据已经用附件上传，经过我个人测试是可以正确下载的。

再次感谢您的帮忙。

fcode

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

fcode 发表于 2015-3-25 22:06
Running...
Time =        7.8090000

试验以后难以表达我的喜悦之感!!! 现在我可以用fs为单位进行计算了(对精确的追求方面有些神经质)。您说得对，这个是correlation function (绝大多数专业词我不会用中文表达，非常抱歉)相关的计算。就是Satchler那本Anguler Momentum147页提到的closure部分的验证。这个角度cos(theta)是stochastic quantity里面的一个时间相关的参数(我表达的不好，可能词不达意)。为了统计上的精确，不但时间间隔取1，还要取2，取3.....假设时间为10，那么要取到至少间隔5; 就是用第一个数据和第六个数据计算，第二个数据和第七个数据计算.....以此类推。就结果本身来说其实没有特别大的差别，但如果plot出来就会发现如果不这样图像会出现锯齿，而不是smooth的。
---
我能不能请教您，我其实看不太出来究竟是哪里的修改起了关健性的作用。为什么计算时间会相差上百倍? 能不能请您简单地谈一谈，根据经验，如何能在做这类步骤很多的重复的计算时优化速度。非常感谢!!!

fcode

我想你之前的代码可能最严重的错误就是
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 在中文里叫自相关。
如果是一维的情况，数据量大一些我会选择在频率域进行。
不过是否适合你的三维情况，还需要推导一下才知道。