希萊雅桑 发表于 2017-4-3 01:42:32

關於超出邊界的問題

程式檔案txt檔:(我用.f跑)

程式內容:
               program sandpile
                integer Z,h(3000,3000),L,t,ix,iy,an(3000,3000),a,n,na,aa
                  integer ax1,ay1,ah1(3000,3000),ax2,ay2,ah2(3000,3000)
                  integer ax3,ay3,ah3(3000,3000)
                        integer tn(100000),tx(100000),ty(100000),tx1(100000)
                        integer ty1(1000000),m,i,j,noapi,k,m1,n1,n2,n3,n4
                  integer a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14
                  integer a15,a16,a17,a18
                integer na1,na2,na3,na4,na5,na6,na7,na8,na9,na10,na11
                        integer na12,na13,na14,na15,na16,na17,na18,na19,na20
                        integer na21,na22,na23
                      call CPU_time(time_begin)
                Z=2
                L=32
                  km=8
                  noap=1000000
                Do ix=1,L !設ix平面從一到L
                Do iy=1,L !設iy平面從一到L
                h(ix,iy)=rand()*Z+1 !晶格內的高度1~4隨機
                enddo   
                enddo               
                     p1=0
                     p2=0
                         p3=0
                     p4=0
                Do noapi=1,noap
      if(km.gt.1)then
                Do ix=1,L
                Do iy=1,L
                  ax1=ix
                  ay1=iy
                  ah1(ax1,ay1)=h(ix,iy)
                enddo   
                enddo
                        mmx=0
      mmx1=0
      na=0
      aa=0
      endif
                  do mm=1,km
                        n=0
                  a=0
                  Do x=1,L
                Do y=1,L
                        an(x,y)=0
                           enddo   
                enddo   
               ix=rand()*L+1
               iy=rand()*L+1
               t=0
               h(ix,iy)=h(ix,iy)+1
                  if(km.gt.1)then            
               if(h(ix,iy).le.Z.and.mmx.eq.1) then
                            kfc=0
                           kfc=rand()*2+1
                   if(kfc.eq.2) mmx=0
      endif
                           if(h(ix,iy).le.Z.and.mmx.eq.0) then
                   do ix=1,L
                   do iy=1,L
      ax2=ix
      ay2=iy
      ah2(ax2,ay2)=h(ix,iy)
      enddo
      enddo
      mmx=1
      endif
      endif
                           if(h(ix,iy).gt.Z) then
                  n=n+1
                        an(ix,iy)=an(ix,iy)+1
                  if(an(ix,iy).eq.1) then
                              a=a+1
                              else
                  a=a+0
                  end if
                  h(ix,iy)=h(ix,iy)-Z
                  t=1
                  tn(1)=1
                  tx(1)=ix
                  ty(1)=iy
101               m=0
                   DO k=1,tn(t)
                      ixtmp=tx(k)
                      iytmp=ty(k)

                   DO 90 kn=1,Z
                  
                        call neighbor(noapi,ixtmp,iytmp,kn,L,i,j,nd)

                     if(nd.eq.1) goto 90

                                 h(i,j)=h(i,j)+1
                     if(h(i,j).gt.Z) then
                         n=n+1
                                       an(i,j)=an(i,j)+1
                     if(an(i,j).eq.1) then
                                 a=a+1
                                 else
                     a=a+0
                     end if
                     m=m+1
                     tx1(m)=i
                      ty1(m)=j
                     h(i,j)=h(i,j)-Z
                     endif
90                   continue
                      enddo
                  if(m.ge.1) then
                      t=t+1
                      tn(t)=m
                      DO m1=1,m
                     tx(m1)=tx1(m1)
                     ty(m1)=ty1(m1)
                      enddo
                        goto 101
                  endif               
                   endif
      if(km.gt.1)then
      if(mmx1.eq.1.and.mmx.eq.0) then
      if(aa.gt.a.and.mmx1.eq.1)mmx1=0
      if(aa.eq.a.and.mmx1.eq.1)then
      kfc1=0
      kfc1=rand()*2+1
      if(kfc1.eq.2) mmx1=0
      endif
      endif
      if(mmx1.eq.0.and.mmx.eq.0) then
                     do ix=1,L
                   do iy=1,L
      ax3=ix
      ay3=iy
      ah3(ax3,ay3)=h(ix,iy)
      enddo
      enddo
      na=n
      aa=a
      mmx1=1
      endif
               Do ax1=1,L
                Do ay1=1,L
                  ix=ax1
                  iy=ay1
                  h(ix,iy)=ah1(ax1,ay1)
                enddo   
                enddo
                        endif      
                     enddo
       if(km.gt.1)then
      if(mmx.eq.1) then
                do ax2=1,L
      do ay2=1,L
      ix=ax2
      iy=ay2
      h(ix,iy)=ah2(ax2,ay2)
      enddo
      enddo
      a=0
      n=0
      endif
      if(mmx1.eq.1.and.mmx.eq.0) then
      do ax3=1,L
      do ay3=1,L
      ix=ax3
      iy=ay3
      h(ix,iy)=ah3(ax3,ay3)
      enddo
      enddo
      a=aa
      n=na
      endif
      endif
                        n1=0
                      n2=0
      n3=0
      n4=0
                      DO i1=1,L
                      DO j1=1,L
                           if(h(i1,j1).eq.1) n1=n1+1
                     if(h(i1,j1).eq.2) n2=n2+1
                         if(h(i1,j1).eq.3) n3=n3+1
                     if(h(i1,j1).eq.4) n4=n4+1
                   enddo
                   enddo
                                       if(2**(0).le.a.and.a.lt.2**(1)) a1=a1+1
                           if(2**(1).le.a.and.a.lt.2**(2)) a2=a2+1
                           if(2**(2).le.a.and.a.lt.2**(3)) a3=a3+1
                           if(2**(3).le.a.and.a.lt.2**(4)) a4=a4+1
                           if(2**(4).le.a.and.a.lt.2**(5)) a5=a5+1
                           if(2**(5).le.a.and.a.lt.2**(6)) a6=a6+1
                           if(2**(6).le.a.and.a.lt.2**(7)) a7=a7+1
                           if(2**(7).le.a.and.a.lt.2**(8)) a8=a8+1
                           if(2**(8).le.a.and.a.lt.2**(9)) a9=a9+1
                           if(2**(9).le.a.and.a.lt.2**(10)) a10=a10+1
                           if(2**(10).le.a.and.a.lt.2**(11)) a11=a11+1
                           if(2**(11).le.a.and.a.lt.2**(12)) a12=a12+1
                           if(2**(12).le.a.and.a.lt.2**(13)) a13=a13+1                     
                                       if(2**(13).le.a.and.a.lt.2**(14)) a14=a14+1
                         if(2**(14).le.a.and.a.lt.2**(15)) a15=a15+1
                         if(2**(15).le.a.and.a.lt.2**(16)) a16=a16+1
                         if(2**(16).le.a.and.a.lt.2**(17)) a17=a17+1
                         if(2**(17).le.a.and.a.lt.2**(18)) a18=a18+1
                                       if(2**(0).le.n.and.n.lt.2**(1)) na1=na1+1
                           if(2**(1).le.n.and.n.lt.2**(2)) na2=na2+1
                           if(2**(2).le.n.and.n.lt.2**(3)) na3=na3+1
                           if(2**(3).le.n.and.n.lt.2**(4)) na4=na4+1
                           if(2**(4).le.n.and.n.lt.2**(5)) na5=na5+1
                           if(2**(5).le.n.and.n.lt.2**(6)) na6=na6+1
                           if(2**(6).le.n.and.n.lt.2**(7)) na7=na7+1
                           if(2**(7).le.n.and.n.lt.2**(8)) na8=na8+1
                           if(2**(8).le.n.and.n.lt.2**(9)) na9=na9+1
                           if(2**(9).le.n.and.n.lt.2**(10)) na10=na10+1
                           if(2**(10).le.n.and.n.lt.2**(11)) na11=na11+1
                           if(2**(11).le.n.and.n.lt.2**(12)) na12=na12+1
                           if(2**(12).le.n.and.n.lt.2**(13)) na13=na13+1                     
                                       if(2**(13).le.n.and.n.lt.2**(14)) na14=na14+1
                         if(2**(14).le.n.and.n.lt.2**(15)) na15=na15+1
                         if(2**(15).le.n.and.n.lt.2**(16)) na16=na16+1
                         if(2**(16).le.n.and.n.lt.2**(17)) na17=na17+1
                         if(2**(17).le.n.and.n.lt.2**(18)) na18=na18+1
                           if(2**(18).le.n.and.n.lt.2**(19)) na19=na19+1
                           if(2**(19).le.n.and.n.lt.2**(20)) na20=na20+1
                           if(2**(20).le.n.and.n.lt.2**(21)) na21=na21+1
                           if(2**(21).le.n.and.n.lt.2**(22)) na22=na22+1
                           if(2**(22).le.n.and.n.lt.2**(23)) na23=na23+1
                              p1=p1+real(n1)/real(L*L)
                                  p2=p2+real(n2)/real(L*L)
                              p3=p3+real(n3)/real(L*L)
                                  p4=p4+real(n4)/real(L*L)
                  enddo
                         write(*,*) 'n1=',na1,',','a1=',a1
                         write(*,*) 'n2=',na2,',','a2=',a2
                         write(*,*) 'n3=',na3,',','a3=',a3
                         write(*,*) 'n4=',na4,',','a4=',a4
                         write(*,*) 'n5=',na5,',','a5=',a5
                         write(*,*) 'n6=',na6,',','a6=',a6
                         write(*,*) 'n7=',na7,',','a7=',a7
                         write(*,*) 'n8=',na8,',','a8=',a8
                         write(*,*) 'n9=',na9,',','a9=',a9
                         write(*,*) 'n10=',na10,',','a10=',a10
                         write(*,*) 'n11=',na11,',','a11=',a11
                         write(*,*) 'n12=',na12,',','a12=',a12
                         write(*,*) 'n13=',na13,',','a13=',a13
                         write(*,*) 'n14=',na14,',','a14=',a14
                         write(*,*) 'n15=',na15,',','a15=',a15
                         write(*,*) 'n16=',na16,',','a16=',a16
                                       write(*,*) 'n17=',na17,',','a17=',a17
                                       write(*,*) 'n18=',na18,',','a18=',a18
                         write(*,*) 'n19=',na19
                         write(*,*) 'n20=',na20
                         write(*,*) 'n21=',na21
                         write(*,*) 'n22=',na22
                         write(*,*) 'n23=',na23
                              write(*,*) '1的平均',p1/real(noap)
                              write(*,*) '2的平均',p2/real(noap)
                              write(*,*) '3的平均',p3/real(noap)
                              write(*,*) '4的平均',p4/real(noap)
                   call CPU_time(time_end)
               write(*,*) 'time=',time_end-time_begin,'s'
                  read(*,*)
                   stop
                  end


                      subroutine neighbor(noapi,i1,j1,k,L,i,j,nd)
                   k=rand()*4+1
                                  nd=0
                      if(k.eq.1) then
                        i=i1-1
                        if(i.eq.0) nd=1
                        j=j1
                      endif
                  if(k.eq.2) then
                        i=i1
                        j=j1+1
                        if(j.eq.L+1) nd=1
                      endif
                      if(k.eq.3) then
                        i=i1+1
                        if(i.eq.L+1) nd=1
                        j=j1
                      endif
                  if(k.eq.4) then
                        i=i1
                        j=j1-1
                        if(j.eq.0) nd=1
                      endif
                      return
                      end


li913 发表于 2017-4-3 11:00:51

a1,a2 这些变量,没有初始化。

希萊雅桑 发表于 2017-4-3 11:06:41

我是想要他們累加,所以沒初始化,沒初始化會出問題嗎?
我修改看看,謝謝

希萊雅桑 发表于 2017-4-3 11:56:27

數值初始化了,還是有問題耶

Villain 发表于 2017-4-3 16:27:45

楼主可以用 /traceback /check:bounds 来看到底哪个数组出问题了

希萊雅桑 发表于 2017-4-3 16:41:49

請問要加在哪裡?

希萊雅桑 发表于 2017-4-3 22:48:09

本帖最后由 希萊雅桑 于 2017-4-3 22:49 编辑

Villain 发表于 2017-4-3 16:27
楼主可以用 /traceback /check:bounds 来看到底哪个数组出问题了
是在這裡檢查嗎?

希萊雅桑 发表于 2017-4-3 22:50:14

li913 发表于 2017-4-3 11:00
a1,a2 这些变量,没有初始化。

我是想要他們累加,所以沒初始化
初始化了還是有問題

希萊雅桑 发表于 2017-4-3 22:52:18

Villain 发表于 2017-4-3 16:27
楼主可以用 /traceback /check:bounds 来看到底哪个数组出问题了

我的km設在4沒有這個問題,設在8就有問題了,不知道是哪裡出界了,可是檢查程式碼都有正確壓QAQ:'(:-dizzy:

Villain 发表于 2017-4-4 09:12:11

希萊雅桑 发表于 2017-4-3 22:52
我的km設在4沒有這個問題,設在8就有問題了,不知道是哪裡出界了,可是檢查程式碼都有正確壓QAQ:-dizz ...

看了看你的代码,貌似你是在用随机函数。debug了一下,程序报错的位置确实是rand所在的行。
就我所知,fortran中没有rand这个随机函数,fortran产生随机数的方法是call random_number(x),为了让随机数更均匀,一般用 call random_seed()来提供种子。
页: [1] 2
查看完整版本: 關於超出邊界的問題