fortran程序报错讨论 - Fortran Coder

[Fortran] 纯文本查看 复制代码
subroutinernde(ua,n)

 c           random number generator proposed bymarsaglia and zaman

 c           in reportfsu-scri-87-50

 c           modified by f. james, 1988and 1989, to generate a vector

 c           of pseudorandom numbers ua oflength n.

 c           modified by k.k

 c    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 c     calling sequences for rnde:                                    ++

 c         call rnde(ua, n)         returns a vector ua of n          ++

 c                      32-bit randomfloating point numbers between  ++

 c                      zero and one.                                 ++

 c         call rnd3i(i1)          initializes the generator from one++

 c                      32-bit integeri1

 c   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

       implicit none

       integer n

       real*8 ua(n), sina(102)

       logical first

       save  first

       integer ijkl, ijklin

       data first/.true./

 c

       if(first) then

          first=.false.

 c           default initialization. userhas called rnde without rnd3i.

          ijkl = 54217137









c     %%%%%%%%%%%%%%%%%%%%%%%

          call rnd3ix(ijkl)



c    %%%%%%%%%%%%%%%%%%%%%%%%%%%









      endif







c     %%%%%%%%%%%%%%%%%%%%%%%%%



      call rnd3x(ua, n)   



c     %%%%%%%%%%%%%%%%%%%%%%%%%%%









      return

 c         initializing routine for rnde,may be called before

 c         generating pseudorandom numberswith rnde. the input

 c         values should be in theranges:  0<=ijklin<=900 oooooo

 c     **************

       entry rnd3i(ijklin)

 c     *************

       first=.false.

       call rnd3ix(ijklin)

       return

 c     ************

       entry rnd3r(sina)

 c     ************

          first=.false.

          call rnd3rx(sina)

       end

       subroutine  rnd3ix(ijkl)

 c             the standard values inmarsaglia's paper, ijkl=54217137

           implicit none

           integer n

           real*8 ua(n), u(97), uni, s, t,zuni

           real*8 sina(102),sout(102)

           integer  jj, m

           integer ns, ijkl, i97, j97, ij,kl, i, j, k, l, ii, ivec

           real*8 twom24, c, cd, cm

           parameter (ns=24,twom24=2.**(-24))

           save c, cd, cm, i97, j97

 c

           ij = ijkl/30082

           kl = ijkl - 30082*ij

           i = mod(ij/177, 177) + 2

           j = mod(ij, 177) + 2

           k = mod(kl/169, 178) + 1

           l = mod(kl, 169)

            do   ii= 1, 97

               s = 0.

               t = .5

                do   jj= 1, ns

                   m = mod(mod(i*j,179)*k,179)

                   i = j

                   j = k

                   k = m

                   l = mod(53*l+1,169)

                   if(mod(l*m,64) .ge. 32)s = s+t

                   t = 0.5*t

                enddo

               u(ii) = s

            enddo

           c = 362436.*twom24

           cd = 7654321.*twom24

           cm = 16777213.*twom24

           i97 = 97

           j97 = 33

           return

 c     ****************

       entry rnd3x(ua, n)

 c     ****************

        do   ivec= 1, n

           uni = u(i97)-u(j97)

           if(uni .lt. 0.)uni=uni+1.

           u(i97) = uni

           i97 = i97-1

           if(i97 .eq. 0) i97=97

           j97 = j97-1

           if(j97 .eq. 0) j97=97

           c = c - cd

           if(c .lt. 0.) c=c+cm

           uni = uni-c

           if(uni .lt. 0.)uni=uni+1.

           ua(ivec) = uni

 c                 replace exact zeros byuniform distr. *2**-24

           if(uni .eq. 0.) then

               zuni = twom24*u(2)

 c               an exact zero here isvery unlikely, but let's be safe.

               if(zuni .eq. 0.) zuni=twom24*twom24

               ua(ivec) = zuni

           endif

        enddo

       return

 c     ****************** to get currentstatus

       entry rnd3s(sout)

 c     ***********

            do   i=1, 97

               sout(i)=u(i)

            enddo

           sout(98)=c

           sout(99)=cd

           sout(100)=cm

           sout(101)=i97

           sout(102)=j97

           return

 c     ****************  to restore the old status

       entry rnd3rx(sina)

            do   i=1, 97

                u(i)=sina(i)

            enddo

           c=sina(98)

           cd=sina(99)

           cm=sina(100)

           i97=sina(101)

           j97=sina(102)

       end