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

24小时热门版块排行榜    

Znn3bq.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 【求助】问一个FORTRAN程序
汕头大学海洋科学接受调剂
51/1返回列表
查看: 736  |  回复: 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的回帖

snoopyzhao

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

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的回帖

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的回帖
查看全部 7 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 22408 312求调剂 +16 门路摸摸 2026-04-14 17/850 2026-04-14 23:59 by Xurambo2014
[考研] 求调剂学校 +11 不会吃肉 2026-04-13 12/600 2026-04-14 19:38 by Art1977
[考研] 08工学 309分求调剂 +12 Yin DY 2026-04-08 12/600 2026-04-14 17:49 by lhj2009
[考研] 297，工科调剂?河南农业大学本科 +9 河南农业大学-能 2026-04-14 9/450 2026-04-14 16:03 by Art1977
[考研] 296求调剂 +8 汪！？！ 2026-04-09 9/450 2026-04-14 15:15 by 逆水乘风
[考研] 求调剂 +16 雪逢冬 2026-04-10 16/800 2026-04-14 14:27 by 逆水乘风
[考研] 305求调剂 +8 玛卡巴卡boom 2026-04-11 8/400 2026-04-14 09:04 by pengliang8036
[考研] 考研求调剂 +12 子木呐 2026-04-12 13/650 2026-04-14 01:19 by 王珺璞
[考研] 求调剂 +12 璃茉一定上岸 2026-04-10 13/650 2026-04-14 00:08 by Equinoxhua
[考研] B区0809 ，数一英一，290 求调剂 +3 泠潍1111 2026-04-12 4/200 2026-04-13 20:35 by 学员JpLReM
[考研] 一志愿0807 数一英一 313 有没有二轮调剂 +12 emokidd 2026-04-08 13/650 2026-04-13 08:32 by lhj2009
[考研] 22408 352分求调剂 +5 努力的夏末 2026-04-09 5/250 2026-04-12 19:17 by wj165256
[考研] 求调剂 +16 张番茄不炒蛋 2026-04-10 17/850 2026-04-12 13:58 by 熬夜成！
[考研] 材料工程日语考生求调剂 +7 0856?调剂 2026-04-10 7/350 2026-04-11 21:33 by 蓝云思雨
[考研] 086003调剂求助 +21 苏弋万 2026-04-09 22/1100 2026-04-11 20:25 by dongdian1
[考研] 284求调剂 +11 archer.. 2026-04-09 12/600 2026-04-11 20:23 by 蓝云思雨
[考研] 283求调剂 +22 那个噜子 2026-04-09 22/1100 2026-04-11 10:41 by 逆水乘风
[考研] 293求调剂 +6 勇远库爱314 2026-04-08 6/300 2026-04-11 10:08 by zhq0425
[考研] 本科211 工科085400 280分求调剂 可跨专业 +11 LZH（等待调剂中 2026-04-10 11/550 2026-04-11 08:39 by zhq0425
[考研] 083200 初试305分 求调剂 暂不考虑跨专业 +15 Claireyyyy 2026-04-09 15/750 2026-04-09 16:11 by zhuimr
信息提示
关闭
请填处理意见
关闭