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

fcode

kif117 发表于 2015-3-24 17:40
您要有好办法就请您用代码介绍介绍，这样我和其他看帖子的人都能尝试运行一下比较哪种更好。我是职业研究 ...

没有人是职业玩 fortran 的。你，我，pasuka 都不是。

楼上的说话一贯就是这样，习惯了就好了。倒是不必太在意了~~

不管是不是职业的，代码里必要的优化还是要做的，是吧？多听听他人的意见，总是能是自己开阔思维。

fcode

kif117 发表于 2015-3-24 22:32
我就是不喜欢上来就摆出教育人的口气说‘你该多看看书’，‘你应该怎样怎样’，‘你怎么样怎么样’，‘你 ...

是的，学生们肯定不愿意听到“回去看看课本”这种回答。

但也就那种就在教科书上明确写着的东西，却有人来问。或者明明是需要整个章节来讲解的知识点，非要让我来给他解释，那我也只能呵呵了......（不指你）

我曾接触过一些人，可能他们就是不喜欢太客气的说话。比如有一个“老师”，就是不允许我叫他老师，并说：“只有你的博硕士导师可以这样称呼”。或者就是不让我称呼他为“您”，而非要我说“你”。

还有一个老外，也很有意思，动辄骂我“stupid”，在我看来很不礼貌。但是他却说：“如果我也犯了这样低级的错误，我很乐意被指责为 stupid。”

尽管如此，上面提到的两位，我还是很尊重他们，并且也发现他们挺可爱的。

最后，确实在编程方面，大家都是业余，这个网站上貌似没有纯学计算机编程的。不过，既然开始写代码了，就把他尽量写完美。

一段优美的代码，不止是可以更准确更快速地得到你的结果，更可以提高你将来再次使用它的频率和效率。对将来你或别人再次阅读它，也会有帮助。（其实说白了，都是方便你自己）

fcode

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

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

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

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

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

fcode

我建议你把 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

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

fcode

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 这些，会更方便你对一个大程序进行模块划分。

fcode

修改为 TRIM(ADJUSTL(testdata)) // '.dat'

rewind(文件通道号) 就可以了，文件通道号需要事先 Open 语句里指定。

Write (a(1), '(3f8.3)') c(2:4) 写法不对，是你先在22楼给出来的。

fcode

我建议你单独开一个帖子来讨论这个问题。

给出你的输入文件（范例，比如里面只有两行），然后你的输出文件应该是怎么样的？重点说明你的需求。表达明确。

我很努力的从你的代码里看出你的意图，但我感觉比较昏乱，而且可能存在逻辑错误。所以，我对你的需求还不了解。