c===========================================================
c     Performs red-black GS relaxation sweep of 2d model 
c     problem and returns norm of running ('instantaneous') 
c     residuals.
c===========================================================
      double precision function relaxc2(u,rhs,nx,ny,h)

         implicit       none

         include       'mgprob.inc'

         integer        nx,              ny
         real*8         u(nx,ny),        rhs(nx,ny)
         real*8         h

         real*8         hm2,      m4hm2,    rresl
         integer        i,        j,
     &                  isw,      jsw,      pass

         relaxc2 =  0.0d0
         hm2     =  1.0d0 / (h * h)
         m4hm2   = -4.0d0 * hm2

         jsw = 1
          do pass = 1 , 2
            isw = jsw
            do j = 2 , ny - 1
               do i = isw + 1 , nx - 1, 2
                  rresl = hm2 * (u(i+1,j) + u(i-1,j) +
     &                           u(i,j+1) + u(i,j-1) -
     &                           4.0d0 * u(i,j)) +
     &                    sigma * u(i,j)**2 -
     &                    rhs(i,j)
                  u(i,j) = u(i,j) - rresl  / 
     &                       (m4hm2 + 2.0d0 * sigma * u(i,j))
                  relaxc2 = relaxc2 + rresl * rresl
               end do
               isw = 3 - isw
            end do
            jsw = 3 - jsw
         end do
         relaxc2 = sqrt(relaxc2 / (nx * ny))

         return

      end

c===========================================================
c    Computes
c               h  h
c      lopu := L  u
c===========================================================
      subroutine lop(lopu,u,nx,ny,h)
 
         implicit       none

         include       'mgprob.inc'
 
         integer        nx,       ny
         real*8         lopu(nx,ny),   u(nx,ny)
         real*8         h
 
         integer        i,       j
         real*8         hm2
 
         hm2 = 1.0d0 / (h * h)
 
         do j = 2 , ny - 1
            do i = 2 , nx - 1
               lopu(i,j) = hm2 * (u(i+1,j) + u(i-1,j) +
     &                            u(i,j+1) + u(i,j-1) - 
     &                            4.0d0*u(i,j)) +
     &                            sigma * u(i,j)**2
            end do
         end do
         call dmlbs(lopu,0.0d0,nx,ny)
 
         return
      end


c===========================================================
c     2--D half-weighted 2:1 restriction with injection
c     of border values.
c===========================================================
      subroutine dmrshw(uf,uc,nxc,nyc)

         implicit       none

         integer        nxc,                nyc
         real*8         uf(2 * nxc - 1, 2 * nyc - 1),
     &                  uc(nxc,nyc)

         integer        ic,       jc,       if,       jf

         do jc = 2 , nyc - 1
            jf = 2 * jc - 1
            do ic = 2 , nxc - 1
               if = 2 * ic - 1
               uc(ic,jc) = 0.5d0 * (uf(if,jf) + 0.25d0 *
     &                     (uf(if+1,jf) + uf(if-1,jf) +
     &                      uf(if,jf+1) + uf(if,jf-1)))
            end do
         end do
         do ic = 1 , nxc
            uc(ic,1  ) = uf(2 * ic - 1,1  )
            uc(ic,nyc) = uf(2 * ic - 1,2 * nyc - 1)
         end do
         do jc = 1 , nyc
            uc(1  ,jc) = uf(1,  2 * jc - 1)
            uc(nxc,jc) = uf(2 * nxc - 1,2 * jc - 1)
         end do

         return

      end


c===========================================================
c     2 : 1 linear interpolation of coarse 2--function 'uc' 
c     to fine function 'uf' (prolongation operator).
c===========================================================
      subroutine dmlint(uc,uf,nxc,nyc)

         implicit       none

         integer        nxc,      nyc
         real*8         uc(nxc,nyc),
     &                  uf(2 * nxc - 1,2 * nyc - 1)

         integer        nxf,      nyf
         integer        ic,       jc,       if,       jf

         nxf = 2 * nxc - 1
         nyf = 2 * nyc - 1

         do jc = 1 , nyc
            jf = 2 * jc - 1
            do ic = 1 , nxc
               uf(2 * ic - 1,jf) = uc(ic,jc)
            end do
         end do
         do jf = 1 , nyf , 2
            do if = 2 , nxf - 1 , 2
               uf(if,jf) = 0.5d0 * (uf(if+1,jf) + uf(if-1,jf))
            end do
         end do
         do jf = 2 , nyf - 1 , 2
            do if = 1 , nxf
               uf(if,jf) = 0.5d0 * (uf(if,jf+1) + uf(if,jf-1))
            end do
         end do

         return

      end

c===========================================================
c     Calculates exact solution: 
c
c        u(x,y) = sin(l_x pi x) sin(l_y pi y)
c===========================================================
      subroutine clcu2(u,h,nx,ny)

         implicit       none

         include       'mgprob.inc'

         integer        nx,       ny
         real*8         u(nx,ny)

         real*8         h

         integer        i,        j

         do  j = 1 , ny
            do  i = 1 , nx
                  u(i,j) = sin(pi * lx * (i - 1) * h) *
     &                     sin(pi * ly * (j - 1) * h)
            end do
         end do
         call dmlbs(u,0.0d0,nx,ny)

         return

      end

c===========================================================
c     Calculates right hand side corresponding to above 
c     exact solution.
c===========================================================
      subroutine clcrhs2(rhs,h,nx,ny)

         implicit       none

         include       'mgprob.inc'

         integer        nx,       ny
         real*8         rhs(nx,ny)

         real*8         h

         real*8         uij
         integer        i,        j

         do  j = 1 , ny
            do  i = 1 , nx
               uij =
     &                     sin(pi * lx * (i - 1) * h) *
     &                     sin(pi * ly * (j - 1) * h)
               rhs(i,j) = (-pi**2 * (lx**2 + ly**2) +
     &                     sigma * uij) * uij
            end do
         end do
         call dmlbs(rhs,0.0d0,nx,ny)

         return

      end

c===========================================================
c     Initializes grid structure.
c===========================================================
      subroutine gsinit(lev)
 
         implicit       none
 
         include       'mg.inc'
 
         integer        qalloc 

         integer        lev
         integer        l,        size
 
c-----------------------------------------------------------
c        For each level.
c-----------------------------------------------------------
         do l = 1 , lev
c-----------------------------------------------------------
c           Define mesh dimensions and spacing.
c-----------------------------------------------------------
            if( l .eq. 1 ) then
               nx(1)  = 3
               ny(1)  = 3
               h(1)   = 0.5d0 
            else 
               h(l)   = 0.5d0 * h(l-1)
               nx(l)  = 2 * nx(l-1) - 1
               ny(l)  = 2 * ny(l-1) - 1
            end if
c-----------------------------------------------------------
c           "Allocate" storage for level-l unknown 'u',  
c           right-hand-side 'rhs', relative local truncation
c           error estimate 'tau', and three temporary 
c           arrays 't1', 't2' and 't3'.
c-----------------------------------------------------------
            size    = nx(l) * ny(l)
            u(l)    = qalloc(size)
            rhs(l)  = qalloc(size)
            tau(l)  = qalloc(size)
            t1(l)   = qalloc(size)
            t2(l)   = qalloc(size)
            t3(l)   = qalloc(size)
         end do
 
         return
 
      end

c===========================================================
c     "Memory allocation" routine.
c     
c     'qalloc' returns the current "pointer" (index) into
c     the main storage array 'q' and increments the pointer
c     by size, effectively "allocating" a block of that 
c     size.  If memory is exhausted, the routine prints 
c     an error message and EXITS.
c===========================================================

      integer function qalloc(size)
         implicit       none
      
         include       'mg.inc'

         integer        size 

         logical        first
         data           first / .true. /

         save   
   
         if( first ) then
            qptr  = 1
            first = .false.
         end if
         qalloc = qptr
         if( qptr + size .gt. qsize ) then
            write(0,1000) size, qsize
1000        format(' qalloc: Could not allocate block of ',
     &             'size ',i8,'  Total size: ',i8)
            stop
         end if
         qptr = qptr + size

         return
      end

c===========================================================
c     Loads "border" of 2-D array with scalar.
c===========================================================
      subroutine dmlbs(a,s,d1,d2)

         implicit      none

         integer       d1,         d2
         real*8        a(d1,d2)

         real*8        s

         integer       i1,         i2

         do i2 = 1 , d2
            a(1, i2) = s
            a(d1,i2) = s
         end do

         do i1 = 1 , d1
            a(i1,1)  = s
            a(i1,d2) = s
         end do

         return

      end

c===========================================================
c     Load 2-D array with scalar.
c===========================================================
      subroutine dmls(a,s,d1,d2)
 
         implicit       none
 
         integer        d1,       d2
         real*8         a(d1,d2)
         real*8         s
 
         integer        i,        j
 
         do j = 1 , d2
            do i = 1 , d1
               a(i,j) = s
            end do      
         end do
 
         return
 
      end

c===========================================================
c     Matrix-matrix add.
c===========================================================
      subroutine dmma(a1,a2,a3,d1,d2)
 
         implicit       none
 
         integer        d1,         d2
         real*8         a1(d1,d2),  a2(d1,d2), a3(d1,d2)
 
         integer        i,        j
 
         do j = 1 , d2
            do i = 1 , d1
               a3(i,j) = a1(i,j) + a2(i,j)
            end do
         end do
 
         return
 
      end

c===========================================================
c     Matrix-matrix subtract.
c===========================================================
      subroutine dmms(a1,a2,a3,d1,d2)
 
         implicit       none
 
         integer        d1,         d2
         real*8         a1(d1,d2),  a2(d1,d2), a3(d1,d2)
 
         integer        i,        j
 
         do j = 1 , d2
            do i = 1 , d1
               a3(i,j) = a1(i,j) - a2(i,j)
            end do
         end do
 
         return
 
      end

c===========================================================
c     Matrix-matrix copy.
c===========================================================
      subroutine dmcopy(a1,a2,d1,d2)

         implicit       none

         integer        d1,        d2
         real*8         a1(d1,d2), a2(d1,d2)
 
         integer        i,        j
 
         do j = 1 , d2
            do i = 1 , d1
               a2(i,j) = a1(i,j)
            end do
         end do
 
         return
 
      end

c===========================================================
c     Matrix l2-norm.
c===========================================================
      double precision function dmnrm2(a,d1,d2)

         implicit       none

         integer        d1,       d2
         real*8         a(d1,d2)

         integer        i,        j

         dmnrm2 = 0.0d0

         do j = 1 , d2
            do i = 1 , d1
               dmnrm2 = dmnrm2 + a(i,j) * a(i,j)
            end do
         end do
         dmnrm2 = sqrt(dmnrm2 / (d1 * d2))

         return

      end

c===========================================================
c     Residual tracing routine.   Generates cheesy "2-d"
c     display.
c===========================================================
      subroutine trace(l,string,resnrm)
         implicit      none
 
         character*(*)  string
         integer        l
         real*8         resnrm

         integer        indent
         parameter    ( indent = 2 )
 
         character*128  spaces 
         data           spaces / ' ' /

         write(0,1000)  spaces(1:1 + (l-1) * indent),
     &                  l, string, resnrm
 1000    format(a,'Level ',i2,' ',a,', sweep ',
     &          t50,1p,e8.2)

         return
 
      end

c===========================================================
c     'gft' output routine.
c
c     This routine uses 'gft_out_bbox': bbox = bounding 
c     box which specifies the limits of the grid function's 
c     domain.
c
c         bbox := ( xmin, xmax, ymin, ymax )
c
c     On the SGIs, specify the libraries 
c
c         -lbbhutil -lmfhdf -ldf -ljpeg -lz
c 
c     (in that precise order) at load-time when using 
c     'gft_out_bbox' or any other routines from the  
c     'bbhutil' library.
c===========================================================
      subroutine my_gft_out(l)
         implicit      none

         include      'mg.inc'
         include      'mgprob.inc'

         integer       l

         logical       first
         real*8        bbox(4),     time
         integer       shape(2)

         data          first / .true. /

         save

         if( .not. gft_trace ) return

         if( first ) then
            time = 0.0d0
            bbox(1) = 0.0d0
            bbox(2) = 1.0d0
            bbox(3) = 0.0d0
            bbox(4) = 1.0d0
            first = .false.
         else
            time = time + 1.0d0
         end if

         shape(1) = nx(l)
         shape(2) = nx(l)
         call gft_out_bbox('u',time,shape,2,bbox,q(u(l)))

         return
 
      end