c===========================================================
c     Driver/testing routine for mgv ...
c===========================================================

      implicit        none

      real*8          mgv

      integer         stderr,      stdin,     stdout
      parameter     ( stderr = 0,  stdin = 5, stdout = 6 )

      real*8          s2r8
      real*8          r8arg
      integer         s2i4
      integer         i4arg
      character*2     itoc
      character*128   sarg,         catsqz

      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )

      integer         i4_never 
      parameter     ( i4_never = -2 000 000 000 )

      integer         iargc,  qalloc, qleft, qused
      real*8          mgvcyc, dvnrm2, dvlinf

      integer          p,    pceil
      common  / comp / p,    pceil

c-----------------------------------------------------------------------
c     Working storage ... 
c
c     1 Meg = 2**20  
c-----------------------------------------------------------------------
      integer         maxq
      parameter     ( maxq = 16 * 2**20 )
      real*8          q(maxq)

      integer         maxrank
      parameter     ( maxrank = 3 )
      integer         rank
      integer         shape(maxrank),     size
      real*8          vh(maxrank)

      integer         uf,      uxctf,     rhsf,     duf,
     *                mresf,   xf

      integer         lf,      nf,        ncyc,     icyc,
     *                ncycle,  preswf,    pstswf,   presw, 
     *                pstsw,   p_save
      real*8          resnrm,  e2nrm,     einfnrm,   hf

      logical         ltrace
      parameter     ( ltrace = .true. )

      logical         compute_residual
      parameter     ( compute_residual = .true. )

      if( iargc() .lt. 2 ) go to 100
      rank = i4arg(1,-1)
      if( rank .ne. 1  .and.  rank .ne. 2  .and. rank  .ne. 3 ) then
         write(0,*) 'mg0: Only rank-1, 2 and 3 implemented'
         stop
      end if

      lf   = i4arg(2,-1)
      if( lf .lt. 1 )      go to 100
      nf   = 2 ** lf + 1
      ncyc   = i4arg(3,1)
      preswf = i4arg(4,3)
      pstswf = i4arg(5,3)
      presw  = i4arg(6,3)
      pstsw  = i4arg(7,3)
      if( ltrace ) then
         write(0,5000) rank, nf, rank, ncyc, 
     *                 preswf, pstswf, presw, pstsw
5000     format(' Rank ',i1,'  Problem size: ',i3,'^',i1,
     *          '  ',i3,' (',i2,',',i2,',',i2,',',i2,') cycles'/)
      end if

      hf = 1.0d0 / (nf - 1)

      call ivls(shape,nf,rank)
      call dvls(vh,hf,rank)

      call qreset(maxq)

      size = nf**rank
      uf    = qalloc(size)
      uxctf = qalloc(size)
      rhsf  = qalloc(size)
      duf   = qalloc(size)
      mresf = duf
      xf    = duf

      call clcrhs(q(rhsf),rank,shape,vh)
      call clcu(q(uxctf),rank,shape,vh)
      call dvls(q(uf),0.0d0,size)

      ncycle = 0
      p_save = p
      do icyc = 1 , ncyc
         resnrm = mgv (rank,shape,vh,
     *                 q(uf),q(rhsf),
     *                 q(p_save),maxq-p_save+1,
     *                 preswf,pstswf,ncycle)
         call dvvs(q(uxctf),q(uf),q(duf),size)
         e2nrm = dvnrm2(q(duf),size)
         einfnrm = dvlinf(q(duf),size)
         write(0,5050) ncycle, resnrm, e2nrm, einfnrm
5050     format(' Cycle: ',i2/' |rres|_2:',1pE13.6,
     *          '  |uxct - u|_2:',1pE13.6,
     *          '  |uxct - u|_i:',1pE13.6)
         if( compute_residual ) then
            call cmres(q(mresf),q(uf),q(rhsf),rank,shape,vh)
            e2nrm   = dvnrm2(q(mresf),size)
            einfnrm = dvlinf(q(mresf),size)
            write(0,5075) e2nrm, einfnrm
5075        format(' |res|_2: ',1pE13.6,'  |res|_i:',1pE13.6)
         end if
         write(0,*)
      end do
      if( ltrace ) then
         write(0,*) p_save + qused(), ' words of storage used'
         write(0,*) qleft(), ' words of storage available'
      end if

      stop

100   continue
         write(0,*) 'usage: mg0 <rank> <lf> [<ncyc|1> '//
     *              '<preswf|3> <pstswf|3>] [<presw|3> <pstsw|3>]'
      stop

      end

c===========================================================
c     Front end to clcu1, clcu2 and clcu3 ...
c===========================================================
      subroutine clcu(u,rank,shape,vh)

         implicit     none

         integer      rank
         real*8       u(*),          vh(rank)
         integer      shape(rank)

         if(      rank .eq. 1 ) then
            call clcu1(u,vh(1),shape(1))
         else if( rank .eq. 2 ) then
            call clcu2(u,vh(1),shape(1),shape(2))
         else if( rank .eq. 3 ) then
            call clcu3(u,vh(1),shape(1),shape(2),shape(3))
         else 
            write(0,*) 'clcu: rank = ', rank,' not implemented'
         end if

         return

      end

c===========================================================
c     Front end to clcrhs1, clcrhs2 and clcrhs3 ...
c===========================================================
      subroutine clcrhs(rhs,rank,shape,vh)

         implicit     none

         integer      rank
         real*8       rhs(*),        vh(rank)
         integer      shape(rank)

         if(      rank .eq. 1 ) then
            call clcrhs1(rhs,vh(1),shape(1))
         else if( rank .eq. 2 ) then
            call clcrhs2(rhs,vh(1),shape(1),shape(2))
         else if( rank .eq. 3 ) then
            call clcrhs3(rhs,vh(1),shape(1),shape(2),shape(3))
         else 
            write(0,*) 'clcrhs: rank = ', rank,' not implemented'
         end if

         return

      end

c===========================================================
c     Calculates exact solution u(x)  = sin(l_x pi x) 
c     and corresponding right hand side for laplacian(u) = rhs.
c===========================================================
      subroutine clcu1(u,h,nx)
 
         implicit       none
 
         include       'mgprob.inc'
 
         integer        nx
         real*8         u(nx)
 
         real*8         h
 
         integer        i,        j
 
         do  i = 1 , nx
               u(i) = sin(pi * lx * (i - 1) * h) 
            end do  
         call dvlbs(u,0.0d0,nx)
 
         return
 
      end
 
c===========================================================

      subroutine clcrhs1(rhs,h,nx)
 
         implicit       none
 
         include       'mgprob.inc'
 
         integer        nx
         real*8         rhs(nx)
 
         real*8         h
 
         real*8         ui 
         integer        i
 
         do  i = 1 , nx
            ui = sin(pi * lx * (i - 1) * h)
            rhs(i) = -lx * lx * pi * pi * ui
         end do
         call dvlbs(rhs,0.0d0,nx)
 
         return
 
      end

c===========================================================
c     Calculates exact solution u(x,y) = sin(l_x pi x) sin(l_y pi y)
c     and corresponding right hand side for laplacian(u) = rhs.
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===========================================================
 
      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) = -(lx * lx + ly * ly) 
     *                           * pi * pi * uij
            end do
         end do
         call dmlbs(rhs,0.0d0,nx,ny)
 
         return
 
      end

c===========================================================
c     Calculates exact solution 
c         u(x,y) = sin(l_x pi x) sin(l_y pi y) sin(l_z pi z) 
c     and corresponding right hand side for 
c         laplacian(u) = rhs.
c===========================================================
 
      subroutine clcu3(u,h,nx,ny,nz)
 
         implicit       none
 
         include       'mgprob.inc'
 
         integer        nx,       ny,       nz
         real*8         u(nx,ny,nz)
 
         real*8         h
 
         integer        i,        j,        k
 
         do  k = 1 , nz
            do  j = 1 , ny
               do  i = 1 , nx
                  u(i,j,k) = sin(pi * lx * (i - 1) * h) *
     *                       sin(pi * ly * (j - 1) * h) *
     *                       sin(pi * lz * (k - 1) * h)
               end do  
            end do
         end do
         call d3lbs(u,0.0d0,nx,ny,nz)
 
         return
 
      end
 
c===========================================================
 
      subroutine clcrhs3(rhs,h,nx,ny,nz)
 
         implicit       none
 
         include       'mgprob.inc'
 
         integer        nx,       ny,      nz
         real*8         rhs(nx,ny,nz)
 
         real*8         h
 
         real*8         uijk
         integer        i,        j,       k
 
         do  k = 1 , nz
            do  j = 1 , ny
               do  i = 1 , nx
                  uijk =
     *                       sin(pi * lx * (i - 1) * h) *
     *                       sin(pi * ly * (j - 1) * h) *
     *                       sin(pi * lz * (k - 1) * h)
                  rhs(i,j,k) = -(lx * lx + ly * ly + lz * lz)
     *                           * pi * pi * uijk
               end do
            end do
         end do
         call d3lbs(rhs,0.0d0,nx,ny,nz)
 
         return
 
      end
