Fortran Coder

查看: 35652|回复: 33
打印 上一主题 下一主题

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

[复制链接]

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
跳转到指定楼层
楼主
发表于 2015-3-23 20:19:38 | 只看该作者 回帖奖励 |正序浏览 |阅读模式
最近在算rotation matrix elemetns相关的,其中有一条是closure的check。具体数学公式之类不详细解释,换成通俗语言问题是这样的:

文件A: 有160个数据,分别是x1(i),y1(i),z1(i)

文件B: 有160个数据,分别是x2(i),y2(i),z2(i)

经过简单运算,二文件数据之差是 x(i),y(i),z(i)

现有S=3(COS(X)^2-1)/2的计算公式,其中COS(X)^2是由x(i),y(i),z(i)运算得来。^是平方符号

现在已经写好的程序是这样的:
[Fortran] 纯文本查看 复制代码
program main
    implicit none

    integer i,j
    integer ios
    integer m,n
    character(160) :: filename1,filename2,filename3,filename4,tmp 
    integer id
    integer n1,n2
    real x1,y1,z1
    real x2,y2,z2
    real x(10000),y(10000),z(10000)
    real v1,v2,v3
    real a,b,c
    real up,down,coss(10000),s1(10000),s2(10000),s3(10000),avg
      integer GetFileN,ii,jj
    character(len=10)::status


    print*, 'Running...'
      filename1='c6_'
      filename2='c8_'
      do m=1,1   ! number of files to read
      write(tmp,*)m             
      open(10,file=trim(filename1)//trim(adjustl(tmp))
    1//'.dat',status='old' )
      open(20,file=trim(filename2)//trim(adjustl(tmp))
    1//'.dat',status='old' )
      do i=1,GetFileN(10) 
    read(10,*,iostat=ios) n1,x1,y1,z1   !read c6
    read(20,*,iostat=ios) n2,x2,y2,z2   !read c8
    if (ios /=0) then
    exit
    endif
      x(i)=x2-x1    !c8-c6 x
    y(i)=y2-y1    !c8-c6 y
    z(i)=z2-z1    !c8-c6 z
c    print*,x(i),y(i),z(i)
c      write(30,'(3f12.5)')x(i),y(i),z(i)
    a=0
    b=0
    c=1 ! labframe
    up=(x(i)*a+y(i)*b+z(i)*c)*(x(i)*a+y(i)*b+z(i)*c)
    down=(x(i)*x(i)+y(i)*y(i)+z(i)*z(i))*(a*a+b*b+c*c)
    coss(i)=up/down   ! cos(i)^2
c      print*, coss(i)   
    s1(i)=(3*coss(i)-1)/2  !n=0 d00
    s2(i)=-sqrt(1.5)*sqrt(1-coss(i))*sqrt(coss(i)) ! n=1 d10
    s3(i)=sqrt(0.375)*(1-coss(i)) !n=2 d20
c    print*, s1(i) 
c     print*, s2(i)
c    print*, s3(i)
    enddo
    enddo
    print*, 'Done.'
      close(10)
      close(20)
    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

请忽略n1和n2,不需要他们参与计算。并且对应A和B两种文件其实每种有几万个(所以用了

trim(filename1)//trim(adjustl(tmp))
    1//'.dat')

的写法。其他s2和s3是指另外两个需要计算的公式。

因为最后需要对这几万个求出的s1,s2,s3值求平均,用现在的写法计算起来超!慢(这里特别说明一下,因为需要统计上有效,所以实际上还要取不同的间隔。比如文件A1,A2这种间隔为1的,下面还要计算A1,A3这样间隔为2的......直到达到文件总数的一半为止(比方说总共10个文件,实际要重复计算五次,间隔分别是1,2,3,4,5。我之前用过的一种‘土办法’是把每一个n对应的xyz提取出来单独存成另一个文件,用另一个程序计算: 见下面)。所以现在需要超过一天以上的时间完成全部运算!)。我想请教各位有没有别的办法可以实现更快的计算,谢谢!

---土办法---
[Fortran] 纯文本查看 复制代码
print*, 'Please enter the number of files:'
      read*,nn
          do mm=1,nn
        write(tmp,*)mm         
      open(10,file=trim('S (')//trim(adjustl(tmp))//').dat'
    1,status='old')
c      open(10,file='nozero3534.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(up1(num))
    allocate(below1(num))
    allocate(cos21(num))
    allocate(up2(num))
    allocate(below2(num))
    allocate(cos22(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
    enddo 

      do step=1,NUM/2
      summ=0
do j=1,NUM-step    
      B(j,:)=A(step+j,:)
      x1(j)=a(j,1)
    y1(j)=a(j,2)
    z1(j)=a(j,3)
      up1(j)=z1(j)*z1(j)
    below1(j)=x1(j)*x1(j)+y1(j)*y1(j)+z1(j)*z1(j)
    cos21(j)=up1(j)/below1(j)
      x2(j)=b(j,1)
    y2(j)=b(j,2)
    z2(j)=b(j,3)
      up2(j)=z2(j)*z2(j)
    below2(j)=x2(j)*x2(j)+y2(j)*y2(j)+z2(j)*z2(j)
    cos22(j)=up2(j)/below2(j)
      s(j)=cos21(j)*cos22(j)
    summ=s(j)+summ
    enddo
    avg=summ(j)/(num-step)



分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
35#
 楼主| 发表于 2015-3-30 19:01:12 | 只看该作者
fcode 发表于 2015-3-30 08:43
我建议你单独开一个帖子来讨论这个问题。

给出你的输入文件(范例,比如里面只有两行),然后你的输出文件 ...

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

2015

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1572 元
贡献
676 点

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

34#
发表于 2015-3-30 08:43:59 | 只看该作者
我建议你单独开一个帖子来讨论这个问题。

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

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

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
33#
 楼主| 发表于 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的程序的写法,特别向您请教程序问题的所在!

2015

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1572 元
贡献
676 点

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

32#
发表于 2015-3-29 07:59:00 | 只看该作者
修改为 TRIM(ADJUSTL(testdata)) // '.dat'

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

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

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
31#
 楼主| 发表于 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积木那样‘越玩越想玩’。可见,在有指导的情况下学习是很愉快的! 非常感谢大家诚恳的帮助!

130

帖子

10

主题

0

精华

大师

F 币
617 元
贡献
372 点

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

30#
发表于 2015-3-29 01:38:06 | 只看该作者
接力14-20樓的討論, 請kif117參考
test.f90 (1.98 KB, 下载次数: 4)

718

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
632 元
贡献
323 点

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

29#
发表于 2015-3-27 20:34:15 | 只看该作者
说明你的 c6c8_10.dat 里面没有 10 行呗。
Do J = 1, num  说明至少要读 num 行

46

帖子

8

主题

0

精华

熟手

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

718

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
632 元
贡献
323 点

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

27#
发表于 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
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-11-1 07:17

Powered by Tencent X3.4

© 2013-2024 Tencent

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