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

         implicit       none

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

         real*8         hm2,      mjelm1,   rresl
         integer        i,        j

         relaxl2 = 0.0d0
         hm2     = 1.0d0 / (h * h)
         mjelm1   =  h * h / 4.0d0

         do j = 2 , ny - 1
            do i = 2 , nx - 1
               rresl = hm2 * (u(i+1,j) + u(i-1,j) +
     *                        u(i,j+1) + u(i,j-1) -
     *                        4.0d0 * u(i,j)) -
     *                 rhs(i,j)
               u(i,j) = u(i,j) + mjelm1 *  rresl
               relaxl2 = relaxl2 + rresl * rresl
            end do
         end do
         relaxl2 = sqrt(relaxl2 / (nx * ny))

         return

      end

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

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

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

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

         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)) -
     *                    rhs(i,j)
                  u(i,j) = u(i,j) + mjelm1 *  rresl
                  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     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       'prob.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       'prob.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) = -(lx * lx + ly * ly)
     *                           * pi * pi * uij
            end do
         end do
         call dmlbs(rhs,0.0d0,nx,ny)

         return

      end

c-----------------------------------------------------------
c     Loads "border" of 2--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 matrix 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     Load matrix with random values between -1.0d0 
c     and 1.0d0. SGI-specific.
c-----------------------------------------------------------
      subroutine dmrand(a,d1,d2)
 
         implicit       none
 
         real*8         rand

         integer        d1,       d2
         real*8         a(d1,d2)
 
         integer        i,        j
 
         do j = 1 , d2
            do i = 1 , d1
               a(i,j) = 2.0d0 * rand() - 1.0d0
            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     Computes estimate of spectral radius rho(G) of 
c     residual amplification matrix and convergence rate
c     R = log10(1/rho) whose reciprocal estimates number of 
c     iterations required to take norm(residuals) down by an 
c     order of magnitude.  Returns (0.0, 0.0) on 
c     initial call.
c-----------------------------------------------------------
      subroutine clcrhor(resnrm,rho,r)
         implicit     none
         real*8       resnrm, rho, r

c-----------------------------------------------------------
c        Maintains previous residual norm.
c-----------------------------------------------------------
         real*8       presnrm

         logical      first / .true. /
         save         presnrm,    first

         if( first ) then
            first = .false.
            rho = 0.0d0
            r   = 0.0d0
         else
            if( resnrm .ne. 0.0d0 ) then
               rho = resnrm / presnrm
               r   = log10(1.0d0 / rho)
            else
               rho = 0.0d0
               r   = 0.0d0
            end if
         end if
         presnrm = resnrm

         return
      end