c version of 1994-june-9

c This one is the "flagship" version for the whole pobox family

c (previously it was called pobox.faster.for, now called pobox.symm.for)

c FUNCTION CURRENTLY IN USE:  function 4

c     To change it, make changes to function "fld" around line 302
c     ALSO THEN CHANGE LINE  98   where it prints out number of function in use

cc FUNCTION  USED is determined by the parameter `nfunc' set below
cc  This didn't work!

c WARNING: must reset max_internal in invert to be as large as N

c The normalization for the Action to match  S = (1/2) phi box phi 
c is now known, and it turns out that:
c MUST MULTIPLY BY -2 TO GET CORRECT ACTION
c (this correction applies in  2D only, of course; different one will
c be needed in 4D)

c WARNING: thresh is set for double precision now (1E-9)

c july 17: WARNING: there may be trouble with using "basis" as an
c argument to "invert", when it hasn't been declared in the calling
c program! for now have just deleted it everywhere, since "kernel" is
c not yet implemented anyhow.

c july 16: now using my inversion routine rather than book-routine
c which Alan used.

c This version has `real' everywhere, rather than `float' or 'double prec.'
c For double precision, compile it with -r8
c (notice the fortran function "real" is apparently single precision in the
c same sense as  sin cos etc.)
c
c on rainbow compile with -static as well (this can be crucial!)
c
c Apparently suns can get confused if text surpasses column 72 (!)
c but compiling with  -e  cures it on yukawa => presumably rodan too.
c
c The `parameter' statement, for some compilers MUST be in the standard
c Fortran 77 form:   `parameter ( ... )'

c Is use of `max' AND `n' in `invert' confused?  No; it is just the usual
c fortran nuisance: can't have "dynamical array sizes"

c _______________________________________________________________________
c
c                     SOME IMPORTANT VARIABLES
c
c        n = actual number of points
c        l = incidence (causal) or link matrix, depending on context.
c              uses convention that:  l(j,k)=1  when  j precedes k
c     ndim = dimension of spacetime (only ndim=2 used so far)
c     cinv = was once inverse of green's fcn, but in present version,
c            it will be just work-space passed to the inversion routine,
c            except in subroutines.
c        f = vector of field-values (phi)
c    dummy = matrix to be inverted (also called `du' in some subroutines)
c            later it holds the inverse itself.
c   icheck = same as dummy, used for checking inverse (not used now)
c       fr = value of diagonal of dummy (chosen =0 for now)
c      mex = size of arrays allocated
c     nexp = mean number of points if randomize N too. 
c    thresh= threshold for matrix to be considered singular
c   ________________________________________________________________

      PROGRAM POBOX

c define maximum allowed N (which is twice nexp, though this is wasteful;
c better would be mex = nexp + 10*sqrt(nexp), say.)
      parameter (nexp=256)
      parameter (mex=2*nexp)
cc Next parameter says which function will be used (a stopgap until ALL are
cc     done simultaneously
cc      parameter (nfunc=1)
c next parameter defines spacetime dimension
      parameter (ndim=2)
c next parameter is value of diagonal to be used
      parameter (fr=0.0E0,ifr=0)
      parameter (thresh=1.0E-9)
      integer iseed, n, l, i, j
c      integer k
cc      integer icheck(mex,mex)
      dimension l(mex,mex)
      real  f(mex),c(mex,4),cin(mex,mex),sup(mex),determ,rank
      real dummy(mex,mex)
      real act,pi,s,eps,w,supmin,supmax
c      real tax,sum
      common pi,eps

      pi=2*asin(1.0)

c Read in shrinking parameter
      print *, '   '
      print *, 'Function currently in use is number 5'
      print *,'nipe epsilon'
      read *, eps
      print *, 'epsilon is', eps
c Read in seed and print it
cc	iseed=itime(1)
      print *,'input integer seed (6 digits or less)'
      read *, iseed
      print *,' iseed=',iseed
c........Read in number of points:
      print *,'   input number of points'
      read *, n
      print *,'   initial N=', n
      if (n .gt. mex) then
         print *, 'too many points!'
         call halt
      endif

cc......This could be used to generate a random n with Poisson distribution
c10   n=npt(nexp,iseed)
c     if (n.gt.mex) goto 10
c	print *,' number of sprinkled points =',n

c sprinkle points and evaluate field thereon:

      call sprinkle(mex,n,ndim,c,f,iseed)

cc	write (unit=1,fmt=*)' coordinates :'
cc	write (unit=1,fmt=*)((c(i,j),j=1,ndim),i=1,n)
cc	write (unit=1,fmt=*)' field values :'
cc	write (unit=1,fmt=*)(f(i),i=1,n)

c build causal matrix (called `l' herein):

      call setc(mex,n,c,ndim,l)

cc	write (unit=1,fmt=*)' causal matrix :'
cc	write (unit=1,fmt=*)((l(i,j),j=1,n),i=1,n)

c reorder so that l is upper-triangular (note that only l and f are
c relabeled, not the actual coordinates.

      call relab(mex,n,l,f)

cc 	write (unit=1,fmt=*)' relabelled causal matrix :'
cc	write (unit=1,fmt=*)((l(i,j),j=1,n),i=1,n)
cc	write (unit=1,fmt=*)' relabelled field values :'
cc	write (unit=1,fmt=*)(f(i),i=1,n)

c build link matrix if ndim=4:
      if (ndim.eq.4) call link(mex,l,n)

c setting dummy to symmetrized l (except for diagonals, which are set
c to fr respectively ifr)  
      do i=1,n-1
         do j=i+1,n
            dummy(i,j)=l(i,j)
            dummy(j,i)=l(i,j)
         end do
      end do
      do i=1,n
         dummy(i,i)=fr
      end do

c checking for disconnected and non-Hegelian elements
      call discard1(mex,n,dummy,f)
      call discard2(mex,n,dummy,f)
      print *,'  surviving number of points is', n

cc	write (unit=1,fmt=*)' matrix to be inverted :'
cc	write (unit=1,fmt=*)((dummy(i,j),j=1,n),i=1,n)

c....... inverting dummy and putting result in dummy (no longer to cin)
      call invert(mex, n, thresh, dummy, cin, rank, sup, determ)

c write out determinant
      print *, '  determinant of matrix being inverted: ', determ
c compute condition number and print it
      supmax=0
      supmin=abs(sup(1))
      do i=1,n
         w= abs(sup(i))
         if (w.gt.supmax) supmax=w
         if (w.lt.supmin) supmin=w         
      end do
      print *,'  extreme sups are', supmax, supmin
      print *, '   their ratio (condition number) is', supmax/supmin

c      print *, ' thresh was', thresh

cc  checking inverse matrix (which used to be in cinv):
c      tax=0.0
c      do i=1,n
c         do j=1,n
c            sum=0.0
c            do k=1,n
c               sum=sum+dummy(i,k)*icheck(k,j)
c            end do
c            if( (i.eq.j).and.(abs(sum-1.).gt.tax) )tax=abs(sum-1.)
c            if ((i.ne.j).and.(abs(sum).gt.tax)) tax=abs(sum)
c         end do
c      end do

c...........computing Action:
      s=act(mex,n,dummy,f)
      print *,' Action =',s 

      end
	

c this no doubt is returning a poissson-distributed integer.
      INTEGER FUNCTION NPT(nexp,iseed)
      integer nexp, iseed, n, j
      real x, p
      real ran3, pdv
      x=ran3(iseed)
      n=nexp
      p=0
      j=0
 100  p=p+pdv(j,n)
      if (x.lt.p) goto 200
      j=j+1
      goto 100
 200  npt=j
      end
	
      REAL FUNCTION pdv(i,n)
      integer i, n, j
      real f
      f=0
      do j=1,i
         f=f+alog(real(j))
      end do
      pdv=i*alog(real(n))-f-real(n)
      pdv=exp(pdv)
      end

c This sprinkles points.  In 2D they lie in the diamond with corners at
c the coordinates: (x,t)= (0,-1), (-1,0), (1,0), (0,1) 
c [N.B. order of coordinates is x,t rather than reverse] 
c in other words they lie in the interval between the two points with 
c t = -1, 1 on the t-axis.
c we also call fld to evaluate "phi"

      SUBROUTINE sprinkle(mex,n,ndim,coor,field,iseed)
      integer mex, n, ndim, iseed, i
      real pi,x1,x2,x3,x4,x,t,coor(mex,4),field(mex)
      real ran3, fld
      common pi
      do i=1,n
         x1=ran3(iseed)
         x2=ran3(iseed)

         if (ndim.eq.4) then
            x3=ran3(iseed)
            x4=ran3(iseed)
         end if

         if (ndim.eq.2) then
            coor(i,1)=x1-x2
            coor(i,2)=x1+x2-1
         end if

         if (ndim.eq.4) then
            coor(i,4)=2*x4-1
            coor(i,1)=x1**(1.0D0/3)*coor(i,4)
            coor(i,2)=acos(1-2*x2)
            coor(i,3)=2*pi*x3
         end if

         x=coor(i,1)
         t=coor(i,ndim)
         field(i)=fld(x,t)
      end do
      end

c This function evaluates phi(x,t) 
c there are presently several phi's to choose from.
c It is "shrunk away from the border" in this version, by amount eps
      real FUNCTION fld(x,t)
      real x,t,pi,u,v,u1,v1
      real eps
      common pi,eps

      u1=x+t
      v1=x-t

c......here is the shrinking
      u = u1/(1.0-eps)
      v = v1/(1.0-eps)

c......and here the computation of phi (being zero out of the interval)
      if (u.lt. -1.0 .or. u.gt.1.0  .or. v.lt. -1.0 .or. v.gt.1.0) then
         fld=0.0
         
      else
         
c     Function 1: 
c
         fld=(1+cos(pi*u))*(1+cos(pi*v)) 
         
c     Function 2:
c         fld = (u + v) * (1+cos(pi*u))*(1+cos(pi*v))

c     Function 3:
c         fld = (u - v) * (1+cos(pi*u))*(1+cos(pi*v))

c The next two functions are meant to test lorentz invariance. They both
c     should give zero, and are related by a boost.  Hence, ignoring
c     finite size effects, pobox should do as good a job on one of
c     them as on the other.  Of course, should also test boost-related
c     functions which have nonzero actions.

c     Function 4: 
c         fld =   (1-cos(4*pi*u)) * (1+cos( pi*v))

c     Function 5: 
         fld =   (1-cos(2*pi*u)) * (1-cos(2*pi*v))

      endif

cc      print *, 'u,v, phi are', u, v, fld

      end



c.....here was a previous function Alan had used:
c     fld=(u-1)*(u+1)*(v-1)*(v+1)
c     fld=fld**2


c next subroutine builds causal matrix, l,  from the sprinkled points
c l(j,k) is set to be 1 when j precedes k
c   ***** to be safe should INITIALIZE   l( , ) TO ZERO
c  ****** the 4D part unchecked by RDS *******
      SUBROUTINE setc(mex,n,c,ndim,l)
      integer mex,n, ndim, l,i,j
      dimension l(mex,mex)
      real c(mex,4)
      real dx,dy,dz,dt,ds2
      if (ndim.eq.4) then
         do i=1,n-1
            do j=i+1,n
               dx=c(i,1)*sin(c(i,2))*cos(c(i,3))
               dx=dx-c(j,1)*sin(c(j,2))*cos(c(j,3))
               dy=c(i,1)*sin(c(i,2))*sin(c(i,3))
               dy=dy-c(j,1)*sin(c(j,2))*sin(c(j,3))
               dz=c(i,1)*cos(c(i,2))-c(j,1)*cos(c(j,2))
               dt=c(i,4)-c(j,4)
               ds2=dx**2+dy**2+dz**2-dt**2
               if (ds2.gt.0) goto 100
               if (dt.gt.0) l(j,i)=1
               if (dt.lt.0) l(i,j)=1
 100        end do
         end do
      end if
      if (ndim.eq.2) then
         do i=1,n-1
            do j=i+1,n
               dx=c(i,1)-c(j,1)
               dt=c(i,2)-c(j,2)
               ds2=dx**2-dt**2
               if (ds2.gt.0) goto 200
               if (dt.gt.0) l(j,i)=1
               if (dt.lt.0) l(i,j)=1
               if (dt.eq.0) then
                  print *, 
     +                 'miracle in "setc", identical points sprinkled'
                  call halt
               end if
 200        end do
         end do
      end if
      end

c this relabels l and f so that l(j,k) vanishes unless j precedes k
      SUBROUTINE relab(mex,n,l,f)
      integer mex, n, l, i, j, k, m
      dimension l(mex,mex)
      real x,f(mex)
      do i=n,2,-1
 50      do j=1,i-1
            if (l(i,j).eq.0) goto 100
            do k=1,n
               m=l(i,k)
               l(i,k)=l(j,k)
               l(j,k)=m
            end do
            do k=1,n
               m=l(k,i)
               l(k,i)=l(k,j)
               l(k,j)=m
            end do
            x=f(i)
            f(i)=f(j)
            f(j)=x
c     ***************does next statement do anything?
            goto 50
 100     end do
      end do
      end

cc This seems intended to compute the link matrix; but it appears 
cc to be UNFINISHED.
cc in particular  kh  is never used!
cc if use kh must make it an argument so that its dimension can be "dynamic"
cc WARNING: NOT ALL DECLARATIONS TAKEN CARE OF
c      SUBROUTINE link(mex,k,num)
c      dimension k(mex,mex),kh(mex,mex)
c      do i=1,num
c         do j=1,num
c            do l=1,num
c               if ((k(i,l).eq.1).and.(k(l,j).eq.1)) goto 10
c            end do
c            goto 20
c 10         kh(i,j)=0
c 20      end do
c      end do
c      end


cc the next subroutine deletes totally unrelated elements.
      SUBROUTINE discard1(mex,n,du,f)
      integer mex, n, i, j 
c      integer ich(mex,mex)
      real du(mex,mex),f(mex)
      if (n.eq.1) then
         print *,' matrix of rank 1 and so no deletion possible'
         print *,' exiting discard1 immediately'
         goto 40
      end if
      if (n.lt.1) then
         print *,' matrix of rank < 1 not possible !'
         call halt
      end if
      i=0
 10   i=i+1
      j=0
 20   j=j+1
      if (du(i,j).eq.1) goto 30
      if (j.lt.n) goto 20
      print *,'  removing a totally disconnected element; i=', i
      call del(mex,n,du,i,f)
      i=i-1
 30   if (i.lt.n) goto 10
 40   end

cc next subroutine deletes each "non-Hegelian element" it finds; until
cc none remains. (One might first find all in a group and replace them
cc by their barycenter (eh))
      SUBROUTINE discard2(mex,n,du,f)
      integer mex, n,  i, j, k
c      integer ich(mex,mex)
      real du(mex,mex),f(mex)
      if (n.eq.1) then
         print *,' matrix of rank 1 and so no deletion possible'
         print *,' exiting discard 2 immediately'
         goto 50
      end if
      if (n.lt.1) then
         print *,' matrix of rank < 1 not possible !'
         call halt
      end if
      i=0
 10   i=i+1
      j=i
 20   j=j+1
      k=0
 30   k=k+1
      if (du(i,k).ne.du(j,k)) goto 40
      if (k.lt.n) goto 30
      print *,' deleting a "non-Hegelian element"; j =', j
      print *, ' the @ pair of field values:', f(i), f(j)
      call del(mex,n,du,j,f)
      j=j-1
 40   if (j.lt.n) goto 20
      if (i.lt.n-1) goto 10
 50   end


c this subroutine deletes the i-th element, specifically it modifies
c the variables n,  du(mmy) and f.
c It does NOT delete the coordinates themselves.   (nor icheck now)
      SUBROUTINE del(mex,n,du,i,f)
      integer mex, n, i, k, m
c      integer ich(mex,mex)
      real du(mex,mex),f(mex)
c.......check for erroneous arguments:
      if (n.le.1) then
         print *,' matrix of rank < 2 and so will not delete'
         call halt
      end if
      if ((i.lt.1).or.(i.gt.n)) then
         print *,' trying to delete point which does not exist !'
         call halt
      end if

      if (i.eq.n) goto 10

c.......modify du and f:
      do k=i+1,n
         do m=1,n
            du(k-1,m)=du(k,m)
c            ich(k-1,m)=ich(k,m)
         end do
         f(k-1)=f(k)
      end do
c.......now taking care of the columns:
      do m=1,n-1
         do k=i+1,n
            du(m,k-1)=du(m,k)
c            ich(m,k-1)=ich(m,k)
         end do
      end do
 10   n=n-1
      end

c this function computes the Action given the inverse matrix and field values
c **** NORMALIZATION NOT CORRECT, THOUGH.
      REAL FUNCTION act(mex,n,cinv,f)
      integer mex,n,i,j
      real cinv(mex,mex),f(mex),pi,sum
      common pi

      sum=0
      do i=1,n
         do j=1,n
            sum=sum+f(i)*cinv(i,j)*f(j)
         end do
      end do
      act=sum
      end
	

c 1993-mar-25
c THIS IS THE SUBROUTINE INVERT
c ( a fortran modification of the pascal program, invert.pas)
c modified from initial version hoping to speed up by about 50% by
c rearranging columns too, and using that certain elements are known to
c have been made 0 already;  or (for B) are still 0.
c
c Idea is that  MAP=1 implies PMA=1  whence A-inv = PM.
c intent is to minimize O(n^3) stuff, which permuting rows and columns
c is not.
c Basically should give identical output with small exception that 
c can find different pivot elements, in those (unusual) cases where
c two potential pivot elements are exactly of same magnitude.

c  This is a subroutine to invert a matrix, A, and find its determinant,
c  or to find its kernel if it is not invertible.                     
c (Also one can  f o r c e  a kernel of desired dimension.)       
c
c  If A is singular enough (as determined by the "constant"
c  thresh) it will return the rank of A as well as a basis for its
c  kernel.
c  If thresh is a negative integer then it will force a kernel of that
c  dimension (or greater of course)                                    
c
c  Notes on internal operation:
c  Ultimately A is replaced by its inverse, the original being saved
c  in A_0 (which for now is sacrificed to save memory).
c  During most of time the inverse is being computed in B.
c  At step mm, column number mm is reduced to canonical form.
c  After finishing the mm-loop, A holds the identity matrix and B holds
c  a matrix M such that  MAP=1 (A being the original matrix), where P is a 
c  permutation matrix; hence P-inv * M is the desired inverse.
c  If A is non-invertible then the columns of the array `basis' should
c  comprise a basis for Ker(A) [so far not implemented]
c************************************************************************
c Meaning of some variables:
c          max_internal : must reset so that max doesn't exceed it
c	   max = dimension of array holding matrix passed to `invert'
c	     n = actual dimension of matrix to be inverted
c                 (unfortunately need both these parameters to overcome
c                 stupidity of fortran.)
c	thresh = threshold for sup(mm)/sup(1) to be called zero;
c                 if it is < 0 then it enforces a kernel of at least
c                 that dimension (once this is implemented!).
c            A = original matrix, converted to its inverse
c            B = working matrix which becomes permuted inverse (see above).
c       determ = det(A)
c       parity    keeps track of all exchanges (row and column)
c         rank =  the rank of A, if it is singular.
c        basis    holds a basis for ker(A), if A is singular (once 
c                 implemented.)       
c	    rp = vector keeping track of the row permutaions
c           cp = vector keeping track of the column permutations 
c                  (columns of B not permuted, only of A)
c                  (cp is same, essentially, as perm of initial version)
c        point   is a pointer which points to the "non-canonized" columns
c                (I think)(used by kernel, when implemented)
      SUBROUTINE invert(max, n, thresh, A, B, rank, sup, determ)
      parameter (max_internal=1024)
      integer  i, j, jw, i_0, j_0, n, parity, rank, mm, max
c       integer k, point(max)
      integer rp(max_internal), cp(max_internal)
      real thresh, A(max,max), B(max,max), sup(max), determ
c       real A_0(max,max)
      real  w, logdet
      integer sign

c test if dimension exceeds internal-maximum or max
      if ((n .gt. max_internal).or.(n.gt.max))then
         print *, 'dimension of',n,
     +   ' is too great for internal parameter or'
         print *,' declared maximum matrix size ! (which are these:)'
         print *, max_internal, max
         stop
      end if

cc can use this to save original A and then restore it at end
cc set A_0 to A
c      do j=1,n   
c         do i=1,n
c            A_0(i,j)=A(i,j)
c         end do
c      end do

c initialize B to identity matrix (first part may not be needed)
      do j=1,n
         do i=1,n
            B(i,j)=0.0
         end do
      end do
      do i=1,n
         B(i,i) = 1.0
      end do

c initialize parity
      parity= 1

c initialize rp and cp
      do i=1,n
         rp(i) = i
      end do
      do j=1,n
         cp(j) = j
      end do

c begin the reduction:
      do mm=1,n

c     find pivot element, sup[mm], at i_0, j_0
         w = -1
         do i= mm, n
            do  j = mm, n 
               if (w .lt. abs(A(i,j))) then
                  i_0 = i
                  j_0 = j
                  w = abs(A(i,j))
               end if
            end do
         end do 
         sup(mm) = A(i_0,j_0)
         
c check whether only zeroes remain; if so invoke Kernel
         w= abs( sup(mm) )
         if ((w.eq.0)  .or.  
     +        (w .lt. thresh*abs(sup(1)) ) .or.  
     +        (mm.eq.n+thresh+1))  then
            rank = mm-1 
            print *,' matrix is singular with rank', rank,
     +           '  and nullity', n-rank
            print *,' sup is ' 
            do i=1 , n
               print *, i,':   ', sup(i)
            end do 
ccc  here must put the instructions which compute kernel
cc   i.e. the equivalent of procedure `kernel'
ccc   call  kernel 
            print *,' kernel not implemented yet, so stopping'
            stop
         end if

c switching rows mm and i_0 of A and B (and rp) and adjusting parity
         if (i_0 .ne. mm) then

            parity = - parity

            i = rp(mm)
            rp(mm) = rp(i_0)
            rp(i_0) = i
           
            do j = 1, n
               w = A(mm,j)   
               A(mm,j) = A(i_0,j)   
               A(i_0,j) = w  
            end do

            do j = 1 , n 
               w = B(mm,j)   
               B(mm,j) = B(i_0,j)  
               B(i_0,j) = w  
            end do
         end if

c switching columns mm and j_0 of B (and cp) and adjusting parity
         if (j_0 .ne. mm) then

            parity = - parity

            j = cp(mm)
            cp(mm) = cp(j_0)
            cp(j_0) = j
           
            do i = 1, n 
               w = A(i,mm)
               A(i,mm) = A(i,j_0)   
               A(i,j_0) = w  
            end do
         end if

c normalizing row number mm in A (and doing same to B)
         w = A(mm,mm)

         do j= mm, n   
            A(mm,j) = A(mm,j)/w
         end do

         do j= 1, n
            B(mm,j) = B(mm,j)/w
         end do
         
c reducing rows mm to n of A and B
c is the treatment of B hobbled by the way rp is used?
         do i=1, n  
            if (i .ne. mm) then
               w= A(i,mm)
               
               do j=mm, n 
                  A(i,j) = A(i,j) - w * A(mm,j)
               end do
               
               do j=1, mm
                  jw = rp(j)
                  B(i,jw) = B(i,jw) - w * B(mm,jw)
               end do

            end if
         end do
      end do         
c this ends the mm-loop

c now put the true inverse into A
c     (the use of cp rather than cp-inverse here was a pure guess! but
c     must be correct since cp here corresponds to perm in "slow" version.)
      do i=1, n 
         do j= 1, n
            A(cp(i),j) = B(i,j)
         end do
      end do


cc      print *,' parity, rp, cp, and sup are:'
cc      print *, parity
cc      print *, (rp(i),i=1,n)
cc      print *, (cp(i),i=1,n)
cc      print *, (sup(i),i=1,n)


c find determinant from sup and parity
c also find log of det since det is too big for machine
      determ = parity
      sign =  parity
      logdet = 0
      do i=1,n
         determ = determ * sup(i)
         logdet = logdet + log(abs(sup(i)))
         if (sup(i).lt.0) sign = - sign
      end do 

      print *, ' sign of determ =',sign,'   and log|det| =', logdet
      end
c this ends the whole subroutine


      SUBROUTINE HALT
      print *,' execution of program halted'
      stop
      end


c This is the random number generator we're using for now:
c     The seed must not be larger than mseed=1618033
c
c     This is the true integer version; 	
c     Notice that the argument is a fake: the value returned is
c     independent of it, except the first time it's invoked!
c     The true random variable is not idum at all, but the entire
c     array, "ma".
c     Can initially call with a negative integer, but no need to in
c     the current version (and empirically (on rodan at least) it
c     it seems to function same independently of sign of seed.)
c
CCCCCCCCCCCCCCCCC WARNING CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c     It is important that the variables "iff" and "ma" be STATIC,
c     as well as "inext"  and "inextp" and maybe others. This is
c     presumably what the compiler does on suns; HOWEVER, on rainbow
c     you must compile with the option -static to get this (so
c     presumbaly internal variables are by default "automatic" on rainbow.)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c  Someone had made all vbles beginnig with `m' real for some reason!
c  RDS has undone that:
      REAL FUNCTION RAN3(idum)
      implicit real(a-h,o-z)
      implicit integer*4(i-n)
      parameter (mbig=4000000, mseed=1618033, mz=0, fac=1.0D0/mbig)
      dimension ma(55)
      data iff /0/
c ******* next is crucial on rainbow, unless compile with -static
      static inext, inextp, iff, ma
c......... this part is only the initialilzation of the array ma:
      if ((idum.lt.0).or.(iff.eq.0)) then
         iff=1
         mj=mseed-iabs(idum)
c       check for the seed having been too large:
         if (mj.lt.0) then 
            print *, 'seed too big: mj=', mj
            pause
         endif
         mj=mod(mj,mbig)
         ma(55)=mj
         mk=1
         do i=1,54
            ii=mod(21*i,55)
            ma(ii)=mk
            mk=mj-mk
            if (mk.lt.mz) mk=mk+mbig
            mj=ma(ii)
         end do
         do k=1,4
            do i=1,55
               ma(i)=ma(i)-ma(1+mod(i+30,55))
               if (ma(i).lt.mz) ma(i)=ma(i)+mbig
            end do
         end do
         inext=0
         inextp=31
         idum=1
      end if
c end of initialization.
      inext=inext+1
      if (inext.eq.56) inext=1
      inextp=inextp+1
      if (inextp.eq.56) inextp=1
      mj=ma(inext)-ma(inextp)
      if (mj.lt.mz) mj=mj+mbig
      ma(inext)=mj
c check mj is in range:
      if((mj.lt.0).or.(mj.gt.mbig))then
         print *,'ran3 is out of range:', mj*fac
         pause
      endif
      ran3=mj*fac
      return
      end

