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

kif117 · 发表于 2015-3-25 23:23:07

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

kif117 · 发表于 2015-3-26 03:35:51

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 · 发表于 2015-3-26 23:11:50

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 · 发表于 2015-3-27 18:20:15

fcode 发表于 2015-3-27 08:30
1. 关于 22 楼的代码，一个非我的专业的，且没有任何注释的代码，我是无法看懂的。最多在语法层次上看懂， ...

首先我必须感谢您的帮助，像这样可以请教具体的问题细节，似乎变得对编程更加感兴趣起来!

   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)的格式。也不明白原因。

献上诚挚的谢意!

kif117 · 发表于 2015-3-27 20:28:23

楚香饭发表于 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‘。能再次请教您吗?

kif117 · 发表于 2015-3-29 04:06:30

chiangtp 发表于 2015-3-29 01:38
接力14-20樓的討論, 請kif117參考

我单纯地汇报一下在我的计算机上运行的时间:

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积木那样‘越玩越想玩’。可见，在有指导的情况下学习是很愉快的! 非常感谢大家诚恳的帮助!

kif117 · 发表于 2015-3-30 03:40:52

fcode 发表于 2015-3-29 07:59
修改为 TRIM(ADJUSTL(testdata)) // '.dat'

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

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的程序的写法，特别向您请教程序问题的所在!

kif117 · 发表于 2015-3-30 19:01:12

fcode 发表于 2015-3-30 08:43
我建议你单独开一个帖子来讨论这个问题。

给出你的输入文件（范例，比如里面只有两行），然后你的输出文件 ...

好的，我另外发了一个新的帖子，麻烦您指导一下。

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