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

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

kif117

kif117 发表于 2015-3-26 03:35
另外，我参考了您数据分类问题求助的帖子，改写了一下提取坐标的程序(之前我是用IF判断ID然后取它的对应 ...

补充: 当我加入批量计算
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'也无法解决问题。

kif117

fcode 发表于 2015-3-25 22:32
我想你之前的代码可能最严重的错误就是
summ=s(j)+summ
这一句了。

另外，我参考了您数据分类问题求助的帖子，改写了一下提取坐标的程序(之前我是用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

fcode 发表于 2015-3-25 22:32
我想你之前的代码可能最严重的错误就是
summ=s(j)+summ
这一句了。

我这里还有一个相关联的很老的程序(可能是我出生以前出版的)。因为不明白什么是‘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

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

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

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 21:23
我建议你把 9300 个数据传上来。
（编辑框里如果有字数限制，可以用附件上传） ...

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

再次感谢您的帮忙。

kif117

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

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

fcode

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