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

24小时热门版块排行榜    

Znn3bq.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 【求助】问一个FORTRAN程序
71/1返回列表
查看: 726  |  回复: 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返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 296求调剂 +11 汪！？！ 2026-04-10 13/650 2026-04-10 23:16 by Ftglcn90
[考研] 求调剂 +9 张番茄不炒蛋 2026-04-10 10/500 2026-04-10 22:07 by 小小虫瓜
[考研] 303求调剂 +7 SereinQ 2026-04-10 8/400 2026-04-10 21:08 by gong120082
[考研] 071000生物学调剂求助 +17 zzzzwww 2026-04-09 20/1000 2026-04-10 15:55 by 求调剂zz
[考研] 环境专硕调剂 +16 会说话的肘子 2026-04-06 16/800 2026-04-10 10:30 by asy1wn
[考研] 268分085602化学工程调剂 +24 月照花林。 2026-04-09 24/1200 2026-04-10 08:09 by Sammy2
[考研] 调剂 +8 只叙离别辞 2026-04-09 10/500 2026-04-09 20:25 by zl8213662
[考研] 085600材料与化工专硕329 求调剂 +24 额cc 2026-04-06 25/1250 2026-04-09 16:01 by wp06
[考研] 265求调剂 +4 风说她早忘了 2026-04-07 4/200 2026-04-09 13:59 by only周
[考研] 0703化学调剂 348分 +14 唉我超真没招了 2026-04-06 15/750 2026-04-08 19:16 by 我减肥1
[考研] 求调剂 +28 111623 2026-04-04 33/1650 2026-04-08 09:24 by 泽润东方
[考研] 计算机408｜在校多次国家级竞赛获奖｜申请调剂 +4 东山大白鹅 2026-04-05 4/200 2026-04-08 00:18 by chongya
[考研] 生物医药调剂｜SCI中科院三区一作+多项科研成果 +8 likangxing 2026-04-07 11/550 2026-04-08 00:02 by lys0704
[考研] 085100建筑学 寻求跨专业调剂 一志愿南大294分 校级省级国家级奖项若干 踏实肯干 +3 1021075758 2026-04-06 4/200 2026-04-07 09:23 by 蓝云思雨
[考研] 材料与化工363求推荐 +11 zh096 2026-04-04 11/550 2026-04-06 19:14 by guanxin1001
[考研] 材料调剂 +5 小刘同学吖吖 2026-04-06 5/250 2026-04-06 18:34 by sherry_1901
[考研] 0817化学工程与技术求调剂，一志愿中海洋319 +14 lv945 2026-04-04 14/700 2026-04-06 10:20 by 蓝云思雨
[考研] 322求调剂 +3 嗯哼哼恒 2026-04-05 3/150 2026-04-05 19:52 by nepu_uu
[考研] 316求调剂 +9 墨辰_Orion926 2026-04-04 9/450 2026-04-04 21:35 by lbsjt
[考研] 求生物学专业调剂-332分 +5 云朵遛弯指南 2026-04-04 5/250 2026-04-04 10:05 by rzh123456
信息提示
关闭
请填处理意见
关闭