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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 【求助】问一个FORTRAN程序
51/1返回列表
查看: 711  |  回复: 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 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 求调剂 +5 yunziaaaaa 2026-03-01 6/300 2026-03-01 23:57 by ccp273206157
[考研] 材料化工调剂 +12 今夏不夏 2026-03-01 13/650 2026-03-01 23:32 by L135790
[考研] 265分求调剂不调专业和学校有行学上就 +6 礼堂丁真258 2026-02-28 8/400 2026-03-01 22:50 by jian_
[考研] 275求调剂 +3 明远求学 2026-03-01 3/150 2026-03-01 22:29 by 刘兵
[考研] 材料类求调剂 +10 wana_kiko 2026-02-28 12/600 2026-03-01 22:10 by 海嵙Y
[考研] 274求调剂 +3 cgyzqwn 2026-03-01 6/300 2026-03-01 21:24 by cgyzqwn
[考研] 299求调剂 +3 Y墨明棋妙Y 2026-02-28 5/250 2026-03-01 21:01 by tangxiaotian
[考研] 306分材料调剂 +4 chuanzhu川烛 2026-03-01 5/250 2026-03-01 19:48 by 无际的草原
[考研] 0856材料求调剂 +11 hyf hyf hyf 2026-02-28 12/600 2026-03-01 18:57 by 18137688336
[考研] 285求调剂 +8 满头大汗的学生 2026-02-28 8/400 2026-03-01 16:47 by caszguilin
[基金申请] 刚录用，没有期刊号，但是在线可看的论文可以放为代表作吗 10+3 arang1 2026-03-01 3/150 2026-03-01 16:43 by babero
[考研] 311求调剂 +6 亭亭亭01 2026-03-01 6/300 2026-03-01 15:41 by 324616
[考研] 307求调剂 +5 wyyyqx 2026-03-01 5/250 2026-03-01 15:21 by Fff-1
[考研] 295复试调剂 +3 简木ChuFront 2026-03-01 3/150 2026-03-01 14:27 by zzxw520th
[考研] 调剂 +3 简木ChuFront 2026-02-28 3/150 2026-03-01 11:46 by 王伟要上岸啊
[论文投稿] 求助coordination chemistry reviews 的写作模板 10+3 ljplijiapeng 2026-02-27 4/200 2026-03-01 09:07 by babero
[考研] 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~~
[考研] 085600材料工程一志愿中科大总分312求调剂 +8 吃宵夜1 2026-02-28 10/500 2026-02-28 20:27 by L135790
[考研] 276求调剂 +3 路lyh123 2026-02-28 4/200 2026-02-28 19:45 by 路lyh123
信息提示
关闭
请填处理意见
关闭