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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 【求助】问一个FORTRAN程序
71/1返回列表
查看: 705  |  回复: 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返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 290求调剂 +8 材料专硕调剂； 2026-02-28 9/450 2026-03-01 12:46 by 闭眼看蓝天
[考研] 306分材料调剂 +3 chuanzhu川烛 2026-03-01 4/200 2026-03-01 12:32 by houyaoxu
[考研] 302材料工程求调剂 +4 Doleres 2026-03-01 5/250 2026-03-01 11:52 by liqiongjy
[考研] 材料学调剂 +7 提神豆沙包 2026-02-28 9/450 2026-03-01 11:50 by 王伟要上岸啊
[考研] 调剂 +3 简木ChuFront 2026-02-28 3/150 2026-03-01 11:46 by 王伟要上岸啊
[考研] 材料类求调剂 +8 wana_kiko 2026-02-28 8/400 2026-03-01 11:44 by 王伟要上岸啊
[考研] 317一志愿华南理工电气工程求调剂 +6 Soliloquy_Q 2026-02-28 11/550 2026-03-01 11:14 by 歌liekkas
[考研] 寻找调剂 +4 LYidhsjabdj 2026-02-28 4/200 2026-03-01 10:56 by sunny81
[硕博家园] 博士自荐 +6 科研狗111 2026-02-26 10/500 2026-03-01 10:02 by 科研狗111
[论文投稿] Optics letters投稿被拒求助 30+3 luckyry 2026-02-26 4/200 2026-03-01 09:06 by babero
[考研] 298求调剂 +5 axyz3 2026-02-28 5/250 2026-03-01 06:45 by 刘兵
[考研] 272求调剂 +4 田智友 2026-02-28 4/200 2026-03-01 06:43 by 刘兵
[考研] 285求调剂 +6 满头大汗的学生 2026-02-28 6/300 2026-03-01 06:29 by Trying]
[基金申请] 面上模板改不了页边距吧？ +5 ieewxg 2026-02-25 6/300 2026-03-01 00:10 by addressing
[考研] 307求调剂 +4 73372112 2026-02-28 6/300 2026-03-01 00:04 by ll247
[考研] 304求调剂 +3 52hz~~ 2026-02-28 5/250 2026-03-01 00:00 by 52hz~~
[考研] 化工专硕348，一志愿985求调剂 +4 弗格个 2026-02-28 6/300 2026-02-28 22:00 by wang_dand
[考博] 26申博 +4 想申博！ 2026-02-26 4/200 2026-02-28 21:37 by limorning
[考研] 264求调剂 +3 巴拉巴拉根556 2026-02-28 3/150 2026-02-28 21:31 by gaoxiaoniuma
[基金申请] 面上可以超过30页吧？ +12 阿拉贡aragon 2026-02-22 13/650 2026-02-26 22:09 by Hahaxia
信息提示
关闭
请填处理意见
关闭