Fortran Coder

楼主: kif117

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

[复制链接]

1948

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1298 元
贡献
547 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2015-3-25 22:32:58 | 显示全部楼层
我想你之前的代码可能最严重的错误就是
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 在中文里叫自相关。
如果是一维的情况,数据量大一些我会选择在频率域进行。
不过是否适合你的三维情况,还需要推导一下才知道。

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
 楼主| 发表于 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

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
 楼主| 发表于 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!

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
 楼主| 发表于 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'也无法解决问题。

1948

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1298 元
贡献
547 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2015-3-27 08:30:10 | 显示全部楼层
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 这些,会更方便你对一个大程序进行模块划分。

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
 楼主| 发表于 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)的格式。也不明白原因。

献上诚挚的谢意!

709

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
596 元
贡献
305 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 2015-3-27 19:58:56 | 显示全部楼层
看注释
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

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
 楼主| 发表于 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‘。能再次请教您吗?

709

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
596 元
贡献
305 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 2015-3-27 20:34:15 | 显示全部楼层
说明你的 c6c8_10.dat 里面没有 10 行呗。
Do J = 1, num  说明至少要读 num 行

130

帖子

10

主题

0

精华

大师

F 币
617 元
贡献
372 点

贡献勋章管理勋章帅哥勋章元老勋章星光勋章规矩勋章

发表于 2015-3-29 01:38:06 | 显示全部楼层
接力14-20樓的討論, 請kif117參考
test.f90 (1.98 KB, 下载次数: 4)
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-29 06:33

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表