      program mg3_demo 
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        'mg_demo' demonstrates the capabilities of a very simple multigrid
c        solver in computing a three dimensional potential field.  This is
c        a simplified multigrid solver in two important respects:
c
c                (1) it solves only a constant coefficient equation,
c                and that only on a uniform cubical grid,
c
c                (2) it solves only a single equation, representing
c                a scalar field rather than a vector field.
c
c        We chose it for its portability and simplicity, and expect that a
c        supercomputer which can run it effectively will also be able to
c        run more complex multigrid programs at least as well.
c
c              Eric Barszcz                        Paul Frederickson
c                                                      RIACS
c        NASA Ames Research Center            NASA Ames Research Center
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      integer lm, nm, nv, nr, nit
      parameter( lm=5, nm=2+2**lm, nv=nm**3 )
      parameter( nr = (8*(nm**3+nm**2+5*nm-23+7*lm))/7 )
      parameter( nit=4 )
c
      real*8 u(nr),v(nv),r(nr),a(0:3),c(0:3)
      real*8 rnm2, rnmu, old2, oldu
      integer ir(lm), mm(lm)
      integer it, n
c
      external timer
      real timer, time0, time1
c
      a(0) = -8.0D0/3.0D0 
      a(1) =  0.0D0 
      a(2) =  1.0D0/6.0D0 
      a(3) =  1.0D0/12.0D0
c
      c(0) =  -3.0D0/8.0D0
      c(1) =  +1.0D0/32.0D0
      c(2) =  -1.0D0/64.0D0
      c(3) =   0.0D0
c
      n  = 2 + 2**lm
      write(6,7)n-2,n-2,n-2,nit
c
      call setup(lm,ir,mm)
      call zero3(u,n,n,n)
      call zran3(v,n,n,n)
c
      time0 = timer( 0.0 )
c
      call resid(u,v,r,n,a)
      call norm2u3(r,n,n,n,rnm2,rnmu,n,n,n)
           old2 = rnm2
           oldu = rnmu
c
      do 20 it=1,nit
        call mg3P(u,v,r,n,a,c,nv,nr,ir,mm,lm)
        call resid(u,v,r,n,a)
 20   continue
c
      time1 = timer( 0.0 )
c
        call norm2u3(r,n,n,n,rnm2,rnmu,n,n,n)
        write(6,'(A40)')
        write(6,'(A40)')' it    L2-norm(R)         unif-norm(R)'
        write(6,6) 0, old2, oldu
        write(6,6) nit, rnm2, rnmu
      write(6,'(A16,F10.2)')'time for solve = ',time1-time0 
      write(6,'(A16,F10.2)')'giving MFLOP = ',
     >      ((31+2+nit*111.14)*(n-2)**3)*1.0D-6 /(time1-time0)
c
 6    format( i4, 2e19.12)
 7    format(/,' Kernel B:  Solving a Poisson problem on a ',i3,' by ',
     > i3,' by ',i3,' grid,',/,' using ',i1,' multigrid iterations.',/)
c
      stop
      end
*
******  mg3P   implements a simple, constant coefficient version of the
******        MG algorithm FAPIN (fast approximate inverse)
*
      subroutine mg3P(u,v,r,n,a,c,nv,nr,ir,mm,lm)
c
      integer n,nv,nr,lm
      real*8 u(nr),v(nv),r(nr)
      real*8 a(0:3),c(0:3)
      integer ir(lm), mm(lm)
c
      integer k, j
c
      do 50 k= lm, 2, -1
      j=k-1
      call rprj3(r(ir(k)),mm(k),r(ir(j)),mm(j))
 50   continue
c
      k = 1
      call zero3(u(ir(k)),mm(k))
      call psinv(r(ir(k)),u(ir(k)),mm(k),c)
      if( lm .eq. 2 )go to 200
      do 100 k = 2, lm-1
      j = k-1
      call zero3(u(ir(k)),mm(k))
      call interp(u(ir(j)),mm(j),u(ir(k)),mm(k))
      call resid(u(ir(k)),r(ir(k)),r(ir(k)),mm(k),a)
      call psinv(r(ir(k)),u(ir(k)),mm(k),c)
 100  continue
 200  continue
      j = lm - 1
      k = lm
      call interp(u(ir(j)),mm(j),u,n)
      call resid(u,v,r,n,a)
      call psinv(r,u,n,c)
c
      return
      end
*
****** psinv applies an approximate inverse as smoother:  u = u + Cr
*
*      This simple implementation costs  27A + 4M per result, where
*      A and M denote the costs of Addition and Multiplication.  
*      By using several two-dimensional buffers one can reduce this
*      cost to  13A + 4M in the general case, or  11A + 3M when the
*      coefficient c(3) is zero.
*
      subroutine psinv(r,u,n,c)
      integer n
      real*8 u(n,n,n),r(n,n,n),c(0:3)
      integer i3, i2, i1
c
      do 600 i3=2,n-1
      do 600 i2=2,n-1
      do 600 i1=2,n-1
 600  u(i1,i2,i3)=u(i1,i2,i3)
     >      +c(0)*( r(i1,  i2,  i3  ) )
     >      +c(1)*( r(i1-1,i2,  i3  ) + r(i1+1,i2,  i3  )
     >           +  r(i1,  i2-1,i3  ) + r(i1,  i2+1,i3  )
     >           +  r(i1,  i2,  i3-1) + r(i1,  i2,  i3+1) )
     >      +c(2)*( r(i1-1,i2-1,i3  ) + r(i1+1,i2-1,i3  )
     >           +  r(i1-1,i2+1,i3  ) + r(i1+1,i2+1,i3  )
     >           +  r(i1,  i2-1,i3-1) + r(i1,  i2+1,i3-1)
     >           +  r(i1,  i2-1,i3+1) + r(i1,  i2+1,i3+1)
     >           +  r(i1-1,i2,  i3-1) + r(i1-1,i2,  i3+1)
     >           +  r(i1+1,i2,  i3-1) + r(i1+1,i2,  i3+1) )
     >      +c(3)*( r(i1-1,i2-1,i3-1) + r(i1+1,i2-1,i3-1)
     >           +  r(i1-1,i2+1,i3-1) + r(i1+1,i2+1,i3-1)
     >           +  r(i1-1,i2-1,i3+1) + r(i1+1,i2-1,i3+1)
     >           +  r(i1-1,i2+1,i3+1) + r(i1+1,i2+1,i3+1) )
c
         call comm3(u,n)
c
      return
      end
*
******* resid computes the residual:  r = v - Au
*
*      This simple implementation costs  27A + 4M per result, where
*      A and M denote the costs of Addition (or Subtraction) and 
*      Multiplication, respectively.  By using several two-dimensional 
*      buffers one can reduce this cost to  13A + 4M in the general 
*      case, or  10A + 3M when the coefficient a(1) is zero.
*
      subroutine resid(u,v,r,n,a)
      integer n
      real*8 u(n,n,n),v(n,n,n),r(n,n,n),a(0:3)
      integer i3, i2, i1
c
      do 600 i3=2,n-1
      do 600 i2=2,n-1
      do 600 i1=2,n-1
 600  r(i1,i2,i3)=v(i1,i2,i3)
     >      -a(0)*( u(i1,  i2,  i3  ) )
     >      -a(1)*( u(i1-1,i2,  i3  ) + u(i1+1,i2,  i3  )
     >                 +  u(i1,  i2-1,i3  ) + u(i1,  i2+1,i3  )
     >                 +  u(i1,  i2,  i3-1) + u(i1,  i2,  i3+1) )
     >      -a(2)*( u(i1-1,i2-1,i3  ) + u(i1+1,i2-1,i3  )
     >                 +  u(i1-1,i2+1,i3  ) + u(i1+1,i2+1,i3  )
     >                 +  u(i1,  i2-1,i3-1) + u(i1,  i2+1,i3-1)
     >                 +  u(i1,  i2-1,i3+1) + u(i1,  i2+1,i3+1)
     >                 +  u(i1-1,i2,  i3-1) + u(i1-1,i2,  i3+1)
     >                 +  u(i1+1,i2,  i3-1) + u(i1+1,i2,  i3+1) )
     >      -a(3)*( u(i1-1,i2-1,i3-1) + u(i1+1,i2-1,i3-1)
     >                 +  u(i1-1,i2+1,i3-1) + u(i1+1,i2+1,i3-1)
     >                 +  u(i1-1,i2-1,i3+1) + u(i1+1,i2-1,i3+1)
     >                 +  u(i1-1,i2+1,i3+1) + u(i1+1,i2+1,i3+1) )
c
         call comm3(r,n)
c
      return
      end
*
****** rprj3 projects onto the next coarser grid, 
****** using a trilinear Finite Element projection:  s = r' = P r
*
*      This simple implementation costs  27A + 4M per result, where
*      A and M denote the costs of Addition and Multiplication.  
*      By using several two-dimensional buffers one can reduce this
*      cost to  23A + 1M.  (Note, however, that the number of results
*      is about 1/8 that of the other three operators.)
*
      subroutine rprj3(r,mk,s,mj)
      integer mk, mj
      real*8 r(mk,mk,mk), s(mj,mj,mj)
      integer j3, j2, j1, i3, i2, i1
c
      do 100 j3=2,mj-1
        i3 = 2*j3-1
      do 100 j2=2,mj-1
        i2 = 2*j2-1
      do 100 j1=2,mj-1
        i1 = 2*j1-1
 100  s(j1,j2,j3) =
     >      0.5D0 * ( r(i1,  i2,  i3  ) )
     >   + 0.25D0 * ( r(i1-1,i2,  i3  ) + r(i1+1,i2,  i3  )
     >             +  r(i1,  i2-1,i3  ) + r(i1,  i2+1,i3  )
     >             +  r(i1,  i2,  i3-1) + r(i1,  i2,  i3+1) )
     >  + 0.125D0 * ( r(i1-1,i2-1,i3  ) + r(i1+1,i2-1,i3  )
     >             +  r(i1-1,i2+1,i3  ) + r(i1+1,i2+1,i3  )
     >             +  r(i1,  i2-1,i3-1) + r(i1,  i2+1,i3-1)
     >             +  r(i1,  i2-1,i3+1) + r(i1,  i2+1,i3+1)
     >             +  r(i1-1,i2,  i3-1) + r(i1-1,i2,  i3+1)
     >             +  r(i1+1,i2,  i3-1) + r(i1+1,i2,  i3+1) )
     > + 0.0625D0 * ( r(i1-1,i2-1,i3-1) + r(i1+1,i2-1,i3-1)
     >             +  r(i1-1,i2+1,i3-1) + r(i1+1,i2+1,i3-1)
     >             +  r(i1-1,i2-1,i3+1) + r(i1+1,i2-1,i3+1)
     >             +  r(i1-1,i2+1,i3+1) + r(i1+1,i2+1,i3+1) )
c
      call comm3(s,mj)
c
      return
      end
*
****** interp adds the trilinear interpolation of the correction
****** from the coarser grid to the current approximation:  u = u + Qu'
*
*      This simple implementation costs  (27/8)A + (7/8)M per result,
*      where A and M denote the costs of Addition and Multiplication.  
*      By using several two-dimensional buffers one can reduce this
*      cost to  (15/8)A + (7/8)M.
*
      subroutine interp(z,m,u,n)
      real*8 z(m,m,m),u(n,n,n)
      integer m, n
      integer i3, i2, i1
c
      do 400 i3=2,m-1
         do 200 i2=2,m-1
            do 100 i1=2,m-1
            u(2*i1-1,2*i2-1,2*i3-1)=u(2*i1-1,2*i2-1,2*i3-1)
     >      +z(i1,i2,i3)
 100        continue
            do 200 i1=2,m-1
            u(2*i1-2,2*i2-1,2*i3-1)=u(2*i1-2,2*i2-1,2*i3-1)
     >      +0.5D0*(z(i1-1,i2,i3)+z(i1,i2,i3))
 200     continue
         do 400 i2=2,m-1
            do 300 i1=2,m-1
            u(2*i1-1,2*i2-2,2*i3-1)=u(2*i1-1,2*i2-2,2*i3-1)
     >      +0.5D0*(z(i1,i2-1,i3)+z(i1,i2,i3))
 300        continue
               do 400 i1=2,m-1
            u(2*i1-2,2*i2-2,2*i3-1)=u(2*i1-2,2*i2-2,2*i3-1)
     >      +0.25D0*(z(i1-1,i2-1,i3)+z(i1-1,i2,i3)
     >            +z(i1,  i2-1,i3)+z(i1,  i2,i3))
 400  continue
c
      do 800 i3=2,m-1
         do 600 i2=2,m-1
            do 500 i1=2,m-1
            u(2*i1-1,2*i2-1,2*i3-2)=u(2*i1-1,2*i2-1,2*i3-2)
     >      +0.5D0*(z(i1,i2,i3-1)+z(i1,i2,i3))
 500        continue
            do 600 i1=2,m-1
            u(2*i1-2,2*i2-1,2*i3-2)=u(2*i1-2,2*i2-1,2*i3-2)
     >      +0.25D0*(z(i1-1,i2,i3-1)+z(i1,i2,i3-1)
     >            +z(i1-1,i2,i3  )+z(i1,i2,i3  ))
 600     continue
         do 800 i2=2,m-1
            do 700 i1=2,m-1
            u(2*i1-1,2*i2-2,2*i3-2)=u(2*i1-1,2*i2-2,2*i3-2)
     >      +0.25D0*(z(i1,i2-1,i3-1)+z(i1,i2,i3-1)
     >            +z(i1,i2-1,i3  )+z(i1,i2,i3  ))
 700        continue
               do 800 i1=2,m-1
            u(2*i1-2,2*i2-2,2*i3-2)=u(2*i1-2,2*i2-2,2*i3-2)
     >      +0.125D0*(z(i1-1,i2-1,i3-1)+z(i1-1,i2,i3-1)
     >             +z(i1,  i2-1,i3-1)+z(i1,  i2,i3-1)
     >             +z(i1-1,i2-1,i3  )+z(i1-1,i2,i3  )
     >             +z(i1,  i2-1,i3  )+z(i1,  i2,i3  ))
 800  continue
c
        call comm3(u,n)
c
      return 
      end
*
****** comm3 communicates on all borders of an array
*
      subroutine comm3(u,n)
      integer n
      real*8 u(n,n,n)
      integer i3, i2, i1
c
      do 100 i3=2,n-1
        do 100 i2=2,n-1
          u(1,i2,i3) = u(n-1,i2,i3)
          u(n,i2,i3) = u( 2, i2,i3)
 100  continue
c
      do 200 i3=2,n-1
        do 200 i1=1,n
          u(i1,1,i3) = u(i1,n-1,i3)
          u(i1,n,i3) = u(i1, 2, i3)
 200  continue
      do 300 i2=1,n
        do 300 i1=1,n
          u(i1,i2,1) = u(i1,i2,n-1)
          u(i1,i2,n) = u(i1,i2, 2 )
 300  continue
c
      return
      end
*
****** norm2u3 evaluates approximations to the L2 norm and the
****** uniform (or L-infinity or Chebyshev) norm, under the
****** assumption that the boundaries are periodic or zero.  Add the
****** boundaries in with half weight (quarter weight on the edges
****** and eigth weight at the corners) for inhomogeneous boundaries.
*
      subroutine norm2u3(r,n1,n2,n3,rnm2,rnmu,nx,ny,nz)
      integer n1, n2, n3, nx, ny, nz
      real*8 rnm2, rnmu, r(n1,n2,n3)
      real*8 s, a
      integer i3, i2, i1
c
      s=0.0D0
      rnmu = 0.0D0
      do 100 i3=2,n3-1
         do 100 i2=2,n2-1
            do 100 i1=2,n1-1
               s=s+r(i1,i2,i3)**2
               a=abs(r(i1,i2,i3))
               if(a.gt.rnmu)rnmu=a
 100  continue
      rnm2=sqrt(s/float((nx-2)*(ny-2)*(nz-2)))
c
      return
      end
*
****** setup   must be called before mg3P
*
      subroutine setup(lm,ir,mm)
      integer lm, ir(lm), mm(lm)
      integer l, k
c
        mm(lm) = 2+2**lm
        do 100 k = lm-1, 1, -1
          mm(k) = 1 + mm(k+1)/2
 100     continue
            l = lm
c
            ir(l)=1
            do 200 k = l-1, 1, -1
            ir(k)=ir(k+1)+mm(k+1)**3
 200        continue
c
      return
      end
*
***
*
      subroutine zero3(z,n)
      integer n
      real*8 z(n,n,n)
      integer i1, i2, i3
c
      do 100 i3=1,n
      do 100 i2=1,n
      do 100 i1=1,n
 100  z(i1,i2,i3)=0.0D0
c
      return
      end
*
****** zran3  loads +1 at ten randomly chosen points,
******        loads -1 at a different ten random points,
******        and zero elsewhere.
*
      subroutine zran3(z,n1,n2,n3)
      integer n1, n2, n3
      real*8 z(n1,n2,n3)
      integer m, i1, i2, i3
      real*8 x, a
      parameter( m = 10,  a = 5.D0 ** 13 )
      real*8 ten( m, 0:1 )
      integer n, i, j1( m, 0:1 ), j2( m, 0:1 ), j3( m, 0:1 )
c
      x = 314159265.D0
c
      do 100 i3=2,n3-1
      do 100 i2=2,n2-1
        call vranl( n1-2, x, a, z(2,i2,i3) )
 100  continue
c
	do 200 i=1,m
	  ten( i, 1 ) = 0.0D0
	  ten( i, 0 ) = 1.0D0
 200	continue
c
      do 300 i3=2,n3-1
      do 300 i2=2,n2-1
      do 300 i1=2,n1-1
	if( z(i1,i2,i3) .gt. ten(1,1) )then
	  ten(1,1) = z(i1,i2,i3) 
	  j1(1,1) = i1
	  j2(1,1) = i2
	  j3(1,1) = i3
	  call bubble( ten, j1, j2, j3, m, 1 )
	endif
	if( z(i1,i2,i3) .lt. ten( 1, 0) )then
	  ten(1,0) = z(i1,i2,i3) 
	  j1(1,0) = i1
	  j2(1,0) = i2
	  j3(1,0) = i3
	  call bubble( ten, j1, j2, j3, m, 0 )
	endif
 300	continue
c 
        write(6,*)' negative charges at'
	write(6,7)(j1(i,0),j2(i,0),j3(i,0),i=1,m)
        write(6,*)' positive charges at'
	write(6,7)(j1(i,1),j2(i,1),j3(i,1),i=1,m)
 7	format(5(' (',i3,2(',',i3),')'))
c
      do 400 i3=1,n3
      do 400 i2=1,n2
      do 400 i1=1,n1
	z(i1,i2,i3) = 0.0D0
 400  continue
c
	do 500 i=1,m
	  z( j1(i,0), j2(i,0), j3(i,0) ) = -1.0D0
	  z( j1(i,1), j2(i,1), j3(i,1) ) = +1.0D0
 500	continue
c
        call comm3(z,n)
c
      return 
      end
*
***** bubble	does a bubble sort in direction dir
*
	subroutine bubble( ten, j1, j2, j3, m, ind )
	real*8 ten( m, 0:1 )
        integer m, ind, j1( m, 0:1 ), j2( m, 0:1 ), j3( m, 0:1 )
	real*8 dir, temp
	integer i, j_temp
c
	if( ind .eq. 1 )then
	  dir = +1.0D0
	else
	  dir = -1.0D0
	endif
c
	do 100 i=1,m-1
	  if( dir*ten(i,ind) .gt. dir*ten(i+1,ind) )then
	    temp = ten( i+1, ind )
	    ten( i+1, ind ) = ten( i, ind )
	    ten( i, ind ) = temp
c
	    j_temp	   = j1( i+1, ind )
	    j1( i+1, ind ) = j1( i,   ind )
	    j1( i,   ind ) = j_temp
c
	    j_temp	   = j2( i+1, ind )
	    j2( i+1, ind ) = j2( i,   ind )
	    j2( i,   ind ) = j_temp
c
	    j_temp	   = j3( i+1, ind )
	    j3( i+1, ind ) = j3( i,   ind )
	    j3( i,   ind ) = j_temp
c
	  else 
	    return
	  endif
c
 100	continue
c
	return
	end
*
***** timer
*
      function timer(x)
      real x, timer
      real secs,tarray(2)
c
c   
c     sun version of timing routine
cccccccccccccccccccccccccccccccccccccc
      secs = etime(tarray)
      timer = secs
cccccccccccccccccccccccccccccccccccccc
c
c       cray version of timing routine
cccccccccccccccccccccccccccccccccccccc
c      timer = second ()
cccccccccccccccccccccccccccccccccccccc
c
c       sgi version of timing routine
cccccccccccccccccccccccccccccccccccccc
c      secs = secnds(0.0)
      timer = secs
cccccccccccccccccccccccccccccccccccccc
      return
      end
*
***** randlc 
*
      FUNCTION RANDLC (X, A)
C
C   This routine returns a uniform pseudorandom double precision number in the
C   range (0, 1) by using the linear congruential generator
C
C   x_{k+1} = a x_k  (mod 2^46)
C
C   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
C   before repeating.  The argument A is the same as 'a' in the above formula,
C   and X is the same as x_0.  A and X must be odd double precision integers
C   in the range (1, 2^46).  The returned value RANDLC is normalized to be
C   between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
C   the new seed x_1, so that subsequent calls to RANDLC using the same
C   arguments will generate a continuous sequence.
C
C   This routine should produce the same results on any computer with at least
C   48 mantissa bits in double precision floating point data.  
C
C   David H. Bailey     October 26, 1990
C
      IMPLICIT REAL*8 (A-H, O-Z)
      IMPLICIT INTEGER (I-N)
      SAVE KS, R23, R46, T23, T46
      DATA KS/0/
C
C   If this is the first call to RANDLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
C   T23 = 2 ^ 23, and T46 = 2 ^ 46.  These are computed in loops, rather than
C   by merely using the ** operator, in order to insure that the results are
C   exact on all systems.  This code assumes that 0.5D0 is represented exactly.
C
      IF (KS .EQ. 0) THEN
        R23 = 1.D0
        R46 = 1.D0
        T23 = 1.D0
        T46 = 1.D0
C
        DO 100 I = 1, 23
          R23 = 0.5D0 * R23
          T23 = 2.D0 * T23
 100    CONTINUE
C
        DO 110 I = 1, 46
          R46 = 0.5D0 * R46
          T46 = 2.D0 * T46
 110    CONTINUE
C
        KS = 1
      ENDIF
C
C   Break A into two parts such that A = 2^23 * A1 + A2 and set X = N.
C
      T1 = R23 * A
      A1 = AINT (T1)
      A2 = A - T23 * A1
C
C   Break X into two parts such that X = 2^23 * X1 + X2, compute
C   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
C   X = 2^23 * Z + A2 * X2  (mod 2^46).
C
      T1 = R23 * X
      X1 = AINT (T1)
      X2 = X - T23 * X1
      T1 = A1 * X2 + A2 * X1
      T2 = AINT (R23 * T1)
      Z = T1 - T23 * T2
      T3 = T23 * Z + A2 * X2
      T4 = AINT (R46 * T3)
      X = T3 - T46 * T4
      RANDLC = R46 * X
C
      RETURN
      END
C
      SUBROUTINE VRANL (N, X, A, Y)
C
C   This routine generates N uniform pseudorandom double precision numbers in
C   the range (0, 1) by using the linear congruential generator
C
C   x_{k+1} = a x_k  (mod 2^46)
C
C   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
C   before repeating.  The argument A is the same as 'a' in the above formula,
C   and X is the same as x_0.  A and X must be odd double precision integers
C   in the range (1, 2^46).  The N results are placed in Y and are normalized
C   to be between 0 and 1.  X is updated to contain the new seed, so that
C   subsequent calls to VRANL using the same arguments will generate a
C   continuous sequence.  If N is zero, only initialization is performed, and
C   the variables X, A and Y are ignored.
C
C   This routine is the standard version designed for scalar or RISC systems.
C   However, it should produce the same results on any single processor
C   computer with at least 48 mantissa bits in double precision floating point
C   data.  
C
C   David H. Bailey     October 26, 1990
C
      IMPLICIT REAL*8 (A-H, O-Z)
      DIMENSION Y(N)
      SAVE KS, R23, R46, T23, T46
      DATA KS/0/
C
C   If this is the first call to VRANL, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
C   T23 = 2 ^ 23, and T46 = 2 ^ 46.  See comments in RANDLC.
C
      IF (KS .EQ. 0 .OR. N .EQ. 0) THEN
        R23 = 1.D0
        R46 = 1.D0
        T23 = 1.D0
        T46 = 1.D0
C
        DO 100 I = 1, 23
          R23 = 0.5D0 * R23
          T23 = 2.D0 * T23
 100    CONTINUE
C
        DO 110 I = 1, 46
          R46 = 0.5D0 * R46
          T46 = 2.D0 * T46
 110    CONTINUE
C
        KS = 1
        IF (N .EQ. 0) RETURN
      ENDIF
C
C   Break A into two parts such that A = 2^23 * A1 + A2 and set X = N.
C
      T1 = R23 * A
      A1 = AINT (T1)
      A2 = A - T23 * A1
C
C   Generate N results.   This loop is not vectorizable.
C
      DO 120 I = 1, N
C
C   Break X into two parts such that X = 2^23 * X1 + X2, compute
C   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
C   X = 2^23 * Z + A2 * X2  (mod 2^46).
C
        T1 = R23 * X
        X1 = AINT (T1)
        X2 = X - T23 * X1
        T1 = A1 * X2 + A2 * X1
        T2 = AINT (R23 * T1)
        Z = T1 - T23 * T2
        T3 = T23 * Z + A2 * X2
        T4 = AINT (R46 * T3)
        X = T3 - T46 * T4
        Y(I) = R46 * X
 120  CONTINUE
C
      RETURN
      END
