应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 【求助】问一个FORTRAN程序
71/1返回列表
查看: 959  |  回复: 6
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

zyj8119

木虫 (著名写手)

[交流] 【求助】问一个FORTRAN程序 已有3人参与

CODE:
subroutine dsrrt(a,xr,xi,n,m,l,b)
        dimension a(n),xr(m),xi(m),b(n)
        if(abs(a(1))+1.0.eq.1.0)then
       l=0
         write(*,5)
         return
        end if
5     format(1x,'err')
      l=1
      k=m
      is=0
        w=1.0
        do 10 i=1,n
10    b(i)=a(i)/a(1)
20    pp=abs(b(k+1))
      if(pp.lt.1.0e-12)then
        xr(k)=0.0
        xi(k)=0.0
        k=k-1
        if(k.eq.1)then
         xr(k)=-b(2)*w/b(1)
         xi(k)=0.0
         return
        end if
        goto 20
        end if
        q=pp**(1.0/k)
        p=q
        w=w*p
        do 30 i=1,k
        b(i+1)=b(i+1)/q
        q=q*p
30    continue
      x=0.0001
        x1=x
        y=0.2
        y1=y
        g=1.0e+37
        dx=1.0
40    u=b(1)
      v=0.0
        do 50 i=1,k
          p=u*x1
          q=v*y1
          pq=(u+v)*(x1+y1)
          u=p-q+b(i+1)
          v=pq-p-q
50    continue
      g1=u*u+v*v
        if(g1.lt.g)goto 105
        if(is.ne.0)goto 80
60    t=t/1.67
      x1=x-t*dx
        y1=y-t*dy
        if(k.ge.50)then
         p=sqrt(x1*x1+y1*y1)
       q=exp(85.0/k)
        if(p.ge.q)goto 60
      end if
        if(t.ge.1.0e-03)goto 40
        if(g.le.1.0e-18)goto 90
65    is=1
      dd=sqrt(dx*dx+dy*dy)
        if(dd.gt.1.0)dd=1.0
        dc=6.28/(k+4.5)
70    c=0.0
80    c=c+dc
      dx=dd*cos(c)
        dy=dd*sin(c)
        x1=x+dx
        y1=y+dy
        if(c.le.6.29)goto 40
        dd=dd/1.67
        if(dd.gt.1.0e-07)goto 70
90    if(abs(y).le.1.0e-06)then
       p=-x
         y=0.0
         q=0.0
        else
         p=-2.0*x
         q=x*x+y*y
         xr(k)=x*w
       xi(k)=-y*w
         k=k-1
        end if
        do 100 i=1,k
         b(i+1)=b(i+1)-b(i)*p
         b(i+2)=b(i+2)-b(i)*p   
100  continue
      xr(k)=x*w
        xi(k)=y*w
        k=k-1
        if(k.eq.1)then
          xr(k)=-b(2)*w/b(1)
          xi(k)=0.0
          return
        end if
      goto 20
105   g=g1
      x=x1
        y=y1
        is=0
        if(g.le.1.0e-22)goto 90
      u1=k*b(1)
        v1=0.0
        do 110 i=2,k
          p=u1*k
          q=v1*y
          pq=(u1+v1)*(x+y)
          u1=p-q+(k-i+1)*b(i)
          v1=pq-p-q
110   continue
      p=u1*u1+v1*v1
        if(p.le.1.0e-20)goto 65
        dx=(u*u1+v*v1)/p
        dy=(u1*v-v1*u)/p
      t=1.0+4.0/k
        goto 60
        end      
         
     
        dimension a(7),xr(6),xi(6),b(7)
        data a/1.0,-5.0,3.0,1.0,-7.0,7.0,-20.0/
        n=7
        m=6
        call dsrrt(a,xr,xi,n,m,l,b)
        if(l.ne.0)then
        do 400 i=1,m
400   continue  
      write(*,500)i,xr(i),xi(i)
      end if
500    format(1x,'x(',i2,1x,')=',e13.6,2x,'j',2x,e13.6)
      end

[ Last edited by 余泽成 on 2010-9-9 at 18:01 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
好好学习，天天向上。
1楼 2010-08-07 22:52:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
老是数组越界，不知道为什么？
赞一下 回复此楼
好好学习，天天向上。
2楼2010-08-07 22:52:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
★
resonant(金币+1):感谢提供：-） 2010-08-07 23:55:03
zyj8119(金币+4):在寻找，你添加的这个程序好像还是数组越界。 2010-08-08 16:44:09
问题在于你的数组 b 下标变成负值了。你跑一下下面的代码看看，加了几个 write，可以看出问题出在哪里，呵呵……至于怎么改，我不了解你的算法，没有办法……
CODE:
        subroutine dsrrt(a,xr,xi,n,m,l,b)
        dimension a(n),xr(m),xi(m),b(n)
        if(abs(a(1))+1.0.eq.1.0)then
       l=0
         write(*,5)
         return
        end if
5     format(1x,'err')
      l=1
      k=m
      is=0
        w=1.0
        do 10 i=1,n
10    b(i)=a(i)/a(1)
20    pp=abs(b(k+1))
      write(*,*) 'i am at line 16', k, pp
      if(pp.lt.1.0e-12)then
        xr(k)=0.0
        xi(k)=0.0
        k=k-1
        write(*,*) 'i am at line 21', k
        if(k.eq.1)then
         xr(k)=-b(2)*w/b(1)
         xi(k)=0.0
         return
        end if
        goto 20
        end if
        q=pp**(1.0/k)
        p=q
        w=w*p
        do 30 i=1,k
        b(i+1)=b(i+1)/q
        q=q*p
30    continue
      x=0.0001
        x1=x
        y=0.2
        y1=y
        g=1.0e+37
        dx=1.0
40    u=b(1)
      v=0.0
        do 50 i=1,k
          p=u*x1
          q=v*y1
          pq=(u+v)*(x1+y1)
          u=p-q+b(i+1)
          v=pq-p-q
50    continue
      g1=u*u+v*v
        if(g1.lt.g)goto 105
        if(is.ne.0)goto 80
60    t=t/1.67
      x1=x-t*dx
        y1=y-t*dy
        if(k.ge.50)then
         p=sqrt(x1*x1+y1*y1)
       q=exp(85.0/k)
        if(p.ge.q)goto 60
      end if
        if(t.ge.1.0e-03)goto 40
        if(g.le.1.0e-18)goto 90
65    is=1
      dd=sqrt(dx*dx+dy*dy)
        if(dd.gt.1.0)dd=1.0
        dc=6.28/(k+4.5)
70    c=0.0
80    c=c+dc
      dx=dd*cos(c)
        dy=dd*sin(c)
        x1=x+dx
        y1=y+dy
        if(c.le.6.29)goto 40
        dd=dd/1.67
        if(dd.gt.1.0e-07)goto 70
90    if(abs(y).le.1.0e-06)then
       p=-x
         y=0.0
         q=0.0
        else
         p=-2.0*x
         q=x*x+y*y
         xr(k)=x*w
       xi(k)=-y*w
         k=k-1
        end if
        write(*,*) 'i am at line 88', k
        do 100 i=1,k
         b(i+1)=b(i+1)-b(i)*p
         b(i+2)=b(i+2)-b(i)*p   
100   continue
      xr(k)=x*w
        xi(k)=y*w
        k=k-1
        if(k.eq.1)then
          xr(k)=-b(2)*w/b(1)
          xi(k)=0.0
          return
        end if
        write(*,*) 'i am at line 101', k
      goto 20
105   g=g1
      x=x1
        y=y1
        is=0
        if(g.le.1.0e-22)goto 90
      u1=k*b(1)
        v1=0.0
        do 110 i=2,k
          p=u1*k
          q=v1*y
          pq=(u1+v1)*(x+y)
          u1=p-q+(k-i+1)*b(i)
          v1=pq-p-q
110   continue
      p=u1*u1+v1*v1
        if(p.le.1.0e-20)goto 65
        dx=(u*u1+v*v1)/p
        dy=(u1*v-v1*u)/p
      t=1.0+4.0/k
        goto 60
        end      
        program test
        dimension a(7),xr(6),xi(6),b(7)
        data a/1.0,-5.0,3.0,1.0,-7.0,7.0,-20.0/
        n=7
        m=6
        call dsrrt(a,xr,xi,n,m,l,b)
        if(l.ne.0)then
        do 400 i=1,m
400   continue  
      write(*,500)i,xr(i),xi(i)
      end if
500    format(1x,'x(',i2,1x,')=',e13.6,2x,'j',2x,e13.6)
      end
赞一下(1人) 回复此楼
3楼2010-08-07 23:31:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

weibangxin
★
nono2009(金币-1):专业区请勿纯表。谢谢。 2010-08-09 06:43:37
zyj8119(金币+1): 2010-09-09 08:29:21
4楼2010-08-08 00:10:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
我说得很清楚啊，我只能找到哪个地方有问题，但我没有能力解决问题，呵呵……
赞一下 回复此楼
5楼2010-08-08 21:25:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
★
resonant(金币+1):en,goto太多的程序很烦人奈...哈哈 2010-08-08 23:31:08
试着跟了一下，发现太困难了，被里面的 goto 搞晕了？这个程序是你写的吗？如果是，想办法重新理一下结构，尽可能少用 goto……
赞一下(1人) 回复此楼
6楼2010-08-08 22:27:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
此贴结贴。
回复此楼
好好学习，天天向上。
7楼2010-09-09 08:29:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zyj8119 的主题更新
71/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 极限拉扯，如何抉择 +6 Isak0320 2026-10-12 8/400 2026-10-12 20:09 by caixiaoye1
[基金申请] 入职后第一次管课题组的专利缴费，对账才发现官费比别人多掏了一倍 +6 13108017953 2026-10-08 6/300 2026-10-12 20:07 by w3w3y3
[有机交流] OpenAI又攻克数学难题了，离AI做有机合成还有多远？ +4 wangyl_123 2026-10-10 8/400 2026-10-12 18:09 by wangyl_123
[教师之家] 课题组会帮你通过考核么？ +3 好运来～9527 2026-10-10 6/300 2026-10-12 14:55 by 简单化xn
[公派出国] 售一区SCI文章T0P，我:8O.551.O54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 22:34 by XQNyDCa1Tk5n
[找工作] 售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 22:26 by XQNyDCa1Tk5n
[公派出国] 售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急 +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 22:26 by XQNyDCa1Tk5n
[考博] 售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 22:15 by XQNyDCa1Tk5n
[考博] 售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急 +3 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 22:11 by XQNyDCa1Tk5n
[教师之家] 售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急 (EPI+-1)(金币-50) +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 22:06 by XQNyDCa1Tk5n
[基金申请] 售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急 +5 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 22:06 by XQNyDCa1Tk5n
[教师之家] 售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急 (EPI+-1)(金币-50) +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 14:25 by VOMexwnDzOFZ
[教师之家] 售SCI文章，我:8O.5.5.1O.54,科目全,可十急 (EPI+-1)(金币-50) +4 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 14:05 by VOMexwnDzOFZ
[考研] 售SCI文章，我:8O.5.5.1O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 13:51 by VOMexwnDzOFZ
[找工作] 售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急 +3 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 10:54 by VOMexwnDzOFZ
[找工作] 售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急 +4 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 10:34 by VOMexwnDzOFZ
[考博] 招收博士生一名 +3 mandolin 2026-10-10 3/150 2026-10-11 10:21 by 不是山谷！
[公派出国] 售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急 +4 qR2AdCWTVtkw 2026-10-09 4/200 2026-10-11 03:14 by M0zrYVSAnxpE
[有机交流] 各位大佬，路线设计求助 +6 zapen 2026-10-08 8/400 2026-10-10 15:59 by zapen
[有机交流] 我应该选哪条合成路线让实验成功率高一些？ +5 jing816y 2026-10-08 9/450 2026-10-10 09:05 by zapen
信息提示
关闭
请填处理意见
关闭