c=======================================================================
c     ROUTINE: mgv
c
c     Solves Poisson equation 
c
c        Laplacian[u] = rhs
c
c     in 1-, 2- and 3-D using fixed  V-cycle LCS 
c     (linear correction scheme) MG.   
c
c     NOTE:  This routine is a real*8 FUNCTION.  Be sure to declare 
c     it as such in the invoking program unit.
c
c     TYPICAL USAGE:  (I) and (O) denote input and ouptut arguments 
c                     respectively
c
c         real*8        mgv
c
c         integer       rank
c         parameter   ( rank = 3)                    ! Problem dimension (I)
c         integer       nx,       ny,       nz
c         parameter   ( nx = 33,  ny = 33,  nz = 33 ) 
c         real*8        u(nx,ny,nz)                  ! Soln (O)
c         real*8        rhs(nx,ny,nz)                ! Rhs  (I)
c         integer       sizeq
c         parameter   ( sizeq = 3 * nx * ny * nz )   ! Size of work array (I)
c         real*8        q(sizeq)                     ! Work array (I)
c         real*8        vh(rank)                     ! Mesh spacings (I)
c         integer       shape(rank)                  ! Mesh dimensions (I)
c         integer       preswp,    pstswp,   ncycle  ! Multi-grid parameters (I)
c        
c
c         vh(1) = h   ! This example assumes mesh spacings in 
c                     ! all coordinate direction = h
c         vh(2) = h
c         vh(3) = h
c         shape(1) = nx
c         shape(2) = ny
c         shape(3) = nz
c         
c         preswp = 3  ! # of pre-CGC sweeps 
c         pstswp = 3  ! # of post-CGC sweeps
c         ncycle = 1  ! # of V-cycles 
c
c         ! Define rhs array
c                 .
c                 .
c                 .
c
c         ! mgv returns norm of residuals after ncycle (preswp,pstswp)
c         ! v-cycles.
c
c         resid = mgv(rank,shape,vh,u,rhs,q,sizeq,preswp,pstswp,ncycle)
c                 .
c                 .
c                 .
c
c
c     IMPLEMENTATION NOTES:  This routine is a front end to mg1v0, 
c     mg2v0 and mg3v0.  It invokes those routines with
c
c        preswf = presw = preswp
c        pstswf = pstsw = pstswp
c=======================================================================
      double precision function mgv
     *   (rank,shape,vh,uf,rhsf,
     *    q,sizeq,
     *    preswp,pstswp,ncycle)

         implicit      none

         real*8        mg1v0,     mg2v0,     mg3v0

         integer       rank,      sizeq
         integer       shape(rank)
         real*8        vh(rank),  uf(*),     rhsf(*),
     *                 q(sizeq)
         integer       preswp,    pstswp,    ncycle

         logical       ltrace
         parameter   ( ltrace = .false. )

         if( ltrace ) then
            write(0,*) 'mgv: Parameter dump'
            write(0,*) 'rank, preswp, pstswp, ncycle'
            write(0,*)  rank, preswp, pstswp, ncycle
            call ivdump(shape,rank,'h',0)
            call dvdump(vh,rank,'h',0)
         end if

         if(      rank .eq. 1 ) then
            mgv = mg1v0(uf,rhsf,shape(1),                  vh(1),
     *                  q,sizeq,
     *                  preswp,pstswp,preswp,pstswp,ncycle)
         else if( rank .eq. 2 ) then
            mgv = mg2v0(uf,rhsf,shape(1),shape(2),         vh(1),
     *                  q,sizeq,
     *                  preswp,pstswp,preswp,pstswp,ncycle)
         else if( rank .eq. 3 ) then
            mgv = mg3v0(uf,rhsf,shape(1),shape(2),shape(3),vh(1),
     *                  q,sizeq,
     *                  preswp,pstswp,preswp,pstswp,ncycle)
         else 
            write(0,*) 'mgv: Invalid rank: ', rank
               mgv = -1.0d0 
         end if
         return

      end

c-----------------------------------------------------------------------
c
c     Front end to mg1v0 which invokes mg1v0  with 
c
c     preswf = presw = preswp
c     pstswf = pstsw = pstswp
c
c-----------------------------------------------------------------------

      double precision function mg1v
     *   (uf,rhsf,nxf,     hf, q,sizeq,
     *    preswp,pstswp,ncycle)

         implicit      none

         real*8        mg1v0

         integer       nxf,                sizeq
         real*8        uf(nxf),            rhsf(nxf),
     *                 q(sizeq)
         real*8        hf
         integer       preswp,    pstswp,  ncycle

         mg1v = mg1v0(uf,rhsf,nxf,    hf,q,sizeq,
     *                preswp,pstswp,preswp,pstswp,ncycle)
         
         return

      end

c-----------------------------------------------------------------------
c
c     Front end to mg2v0 which invokes mg2v0  with 
c
c     preswf = presw = preswp
c     pstswf = pstsw = pstswp
c
c-----------------------------------------------------------------------

      double precision function mg2v
     *   (uf,rhsf,nxf,nyf, hf, q,sizeq,
     *    preswp,pstswp,ncycle)

         implicit      none

         real*8        mg2v0

         integer       nxf,       nyf,     sizeq
         real*8        uf(nxf,nyf),        rhsf(nxf,nyf),
     *                 q(sizeq)
         real*8        hf
         integer       preswp,    pstswp,  ncycle

         mg2v = mg2v0(uf,rhsf,nxf,nyf,hf,q,sizeq,
     *                preswp,pstswp,preswp,pstswp,ncycle)
         
         return

      end

c-----------------------------------------------------------------------
c
c     Front end to mg3v0 which invokes mg3v0  with 
c
c     preswf = presw = preswp
c     pstswf = pstsw = pstswp
c
c-----------------------------------------------------------------------

      double precision function mg3v
     *   (uf,rhsf,nxf,nyf,nzf, hf, q,sizeq,
     *    preswp,pstswp,ncycle)

         implicit      none

         real*8        mg3v0

         integer       nxf,       nyf,     nzf,      sizeq
         real*8        uf(nxf,nyf,nzf),    rhsf(nxf,nyf,nzf),
     *                 q(sizeq)
         real*8        hf
         integer       preswp,    pstswp,  ncycle

         mg3v = mg3v0(uf,rhsf,nxf,nyf,nzf,hf,q,sizeq,
     *                preswp,pstswp,preswp,pstswp,ncycle)
         
         return

      end

c-----------------------------------------------------------------------
c
c     Basic 1d v-cycle multi grid solver for 
c
c     Laplacian[u(x)] = rhs(x)
c
c     This version allows nxf, nyf which are of form
c     2^p * <small number> + 1
c
c     where <small number> .le. maxn0 defined in the parameter 
c     statement below
c
c-----------------------------------------------------------------------

      double precision function mg1v0
     *   (uf,rhsf,nxf,         hf, q, sizeq,
     *    preswf,pstswf,presw,pstsw,ncycle)

         implicit      none

         real*8        relaxc1
         integer       qalloc,    lilog2
         logical       lp2np1

         integer       nxf,                sizeq
         real*8        uf(nxf),            rhsf(nxf), 
     *                 q(sizeq)
         real*8        hf
         integer       preswf,    pstswf,  presw,    pstsw,
     *                 ncycle
c
c        Multi-grid parameters, data structures and pointers 
c
         real*8        epsic
         parameter   ( epsic = 1.0d-10 )

         integer       lmin,      lmax
         parameter   ( lmin = 1,  lmax = 20 )

         integer       nx(lmin:lmax)
         real*8        h(lmin:lmax)

         integer       u(lmin:lmax),       rhs(lmin:lmax),
     *                 mres(lmin:lmax),    du(lmin:lmax)

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

         integer          p,      pceil
         common  / comp / p,      pceil
c
c        For parameter (nxf) checking 
c
         integer       checkfacts

         logical       badargs

         integer       nfact_x,   nx0

         integer       maxn0
         parameter   ( maxn0 = 17 )
c
c        Locals ...
c
         integer       lxf, 
     *                 l,         lf,        lc,       nl,
     *                 sweep

         real*8        resnrm

         logical       ltrace 
         parameter   ( ltrace = .false. )

         save

         mg1v0 = -1.0d0
c
c        Parameter check ...
c

         if( ltrace ) then
            write(0,*) '>>> mg1v0: Parameter dump: '
            write(0,*) nxf, hf
            write(0,*) preswf, pstswf, presw, pstsw, ncycle
         end if

         badargs = .false.
         if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then
            write(*,*) 'mg1v0: Bad nxf value: ', nxf
            badargs = .true.
         end if
         if( badargs) then
            return
         end if

         nl = nfact_x + 1
         lc = 1
         lf = lc + nl - 1

         nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1

         if( ltrace ) then
            write(*,*) 'mg1v0: ', nl, ' levels.'
            write(*,*) 'mg1v0:  Coarse grid ', nx0
         end if
         
         if( ncycle .eq. 0 ) then
c
c           Set up grid structure, allocate storage ...
c
            call qreset(sizeq)
            nx(lf) = nxf
            h(lf)  = hf
            u(lf)    = imposs
            rhs(lf)  = imposs
            mres(lf) = qalloc(nx(lf))
            du(lf)   = mres(lf)
            do 1000 l = lf - 1 , lc , -1
               nx(l) = (nx(l+1) - 1) / 2 + 1
               h(l) = 2.0d0 * h(l+1)

               u(l)    = qalloc(nx(l))
               rhs(l)  = qalloc(nx(l))
               mres(l) = qalloc(nx(l))
               du(l)   = mres(l)
 1000       continue
            if( ltrace ) then
               call ivdump(nx,nl,'nx',0)
               call dvdump(h,nl,'h',0)
               call ivdump(u,nl,'u pointers',0)
               call ivdump(rhs,nl,'rhs pointers',0)
               call ivdump(mres,nl,'mres pointers',0)
               call ivdump(du,nl,'du pointers',0)
            end if
         end if
c
c        Cycle from fine to coarse ...
c
         do l = lf , lc , -1
            if( l .ne. lf  .or.  ncycle .eq. 0 ) then
               if(      l .eq. lc ) then
 100              continue
                     resnrm = relaxc1(q(u(l)),q(rhs(l)),
     *                                nx(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg1v0: pre sweep',l,resnrm)
                        call cres1_trace('mg1v0: pre lc',60,
     &                                   q(mres(l)),q(u(l)),q(rhs(l)),
     &                                   nx(l),h(l))
                     end if
                  if( resnrm .gt. epsic ) go to 100
               else if( l .eq. lf ) then
                  do sweep = 1 , preswf
                     resnrm = relaxc1(uf,rhsf,
     *                                nx(l),h(l))
                     call cmres1(q(mres(l)),uf,rhsf,
     *                           nx(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg1v0: pre sweep',l,resnrm)
                        call cres1_trace('mg1v0: pre lf',60,
     &                                   q(mres(l)),uf,rhsf,
     &                                   nx(l),h(l))
                     end if
                  end do
               else 
                  do sweep = 1 , presw 
                     resnrm = relaxc1(q(u(l)),q(rhs(l)),
     *                                nx(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg1v0: pre sweep',l,resnrm)
                        call cres1_trace('mg1v0: pre inter',60,
     &                                   q(mres(l)),q(u(l)),q(rhs(l)),
     &                                   nx(l),h(l))
                     end if
                  end do
               end if
            end if
c
c           Set up coarse grid problem ...
c
            if(      l .eq. lf ) then
               call cmres1(q(mres(l)),uf,rhsf,nx(l),h(l))
            else if( l .ne. lc ) then
               call cmres1(q(mres(l)),q(u(l)),q(rhs(l)),
     *                    nx(l),h(l))
            end if
            if( l .ne. lc ) then
               call dvrshw(q(mres(l)),q(rhs(l-1)),
     *                     nx(l-1))
               call dvls(q(u(l-1)),0.0d0,
     *                   nx(l-1))
            end if

         end do 
c
c        Cycle from coarse to fine ...
c

         do l = lc + 1 , lf 
c
c           Interpolate correction ...
c
            call dvlint(q(u(l-1)),q(du(l)),nx(l-1))
            if( l .eq. lf ) then
               call dvva(uf,q(du(l)),uf,nx(l))
            else 
               call dvva(q(u(l)),q(du(l)),q(u(l)),nx(l))
            end if

            if( l .eq. lf ) then
               do sweep = 1 , pstswf
                  resnrm = relaxc1(uf,rhsf,nx(l),h(l))
                  if( ltrace ) then
                     call trace_res('mg1v0: post sweep',l,resnrm)
                     call cres1_trace('mg1v0: post lf',60,
     &                                q(mres(l)),uf,rhsf,
     &                                nx(l),h(l))
                  end if
               end do
            else 
               do sweep = 1 , pstsw 
                  resnrm = relaxc1(q(u(l)),q(rhs(l)),
     *                             nx(l),h(l))
                  if( ltrace ) then
                     call trace_res('mg1v0: post sweep',l,resnrm)
                     call cres1_trace('mg1v0: post other',60,
     &                                q(mres(l)),q(u(l)),q(rhs(l)),
     &                                nx(l),h(l))
                  end if
               end do
            end if

         end do

         ncycle = ncycle + 1
         mg1v0 = resnrm

         return

      end


c-----------------------------------------------------------------------
c
c     Basic 2d v-cycle multi grid solver for 
c
c     Laplacian[u(x,y)] = rhs(x,y)
c
c     on rectangular domain. 
c
c     This version allows nxf, nyf which are of form
c     2^p * <small number> + 1
c
c     where <small number> .le. maxn0 defined in the parameter 
c     statement below
c
c-----------------------------------------------------------------------

      double precision function mg2v0
     *   (uf,rhsf,nxf,nyf,     hf, q, sizeq,
     *    preswf,pstswf,presw,pstsw,ncycle)

         implicit      none

         real*8        relaxc2
         integer       qalloc,    lilog2
         logical       lp2np1

         integer       nxf,       nyf,     sizeq
         real*8        uf(nxf,nyf),        rhsf(nxf,nyf), 
     *                 q(sizeq)
         real*8        hf
         integer       preswf,    pstswf,  presw,    pstsw,
     *                 ncycle
c
c        Multi-grid parameters, data structures and pointers 
c
         real*8        epsic
         parameter   ( epsic = 1.0d-10 )

         integer       lmin,      lmax
         parameter   ( lmin = 1,  lmax = 20 )

         integer       nx(lmin:lmax),      ny(lmin:lmax)
         real*8        h(lmin:lmax)

         integer       u(lmin:lmax),       rhs(lmin:lmax),
     *                 mres(lmin:lmax),    du(lmin:lmax)

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

         integer          p,      pceil
         common  / comp / p,      pceil
c
c        For parameter (nxf, nyf) checking 
c
         integer       checkfacts

         logical       badargs

         integer       nfact_x,   nx0
         integer       nfact_y,   ny0

         integer       maxn0
         parameter   ( maxn0 = 17 )
c
c        Locals ...
c
         integer       lxf,       lyf,    
     *                 l,         lf,        lc,       nl,
     *                 sweep

         real*8        resnrm

         logical       ltrace 
         parameter   ( ltrace = .false. )

         save

         mg2v0 = -1.0d0
c
c        Parameter check ...
c

         if( ltrace ) then
            write(0,*) '>>> mg2v0: Parameter dump: '
            write(0,*) nxf, nyf, hf
            write(0,*) preswf, pstswf, presw, pstsw, ncycle
         end if

         badargs = .false.
         if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then
            write(*,*) 'mg2v0: Bad nxf value: ', nxf
            badargs = .true.
         end if
         if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then
            write(*,*) 'mg2v0: Bad nyf value: ', nyf
            badargs = .true.
         end if
         if( badargs) then
            return
         end if

         nl = min(nfact_x,nfact_y) + 1
         lc = 1
         lf = lc + nl - 1

         nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1
         ny0 = (nyf - 1) / (2 ** (nl - 1)) + 1

         if( ltrace ) then
            write(*,*) 'mg2v0: ', nl, ' levels.'
            write(*,*) 'mg2v0:  Coarse grid ', nx0, ' x ', ny0
         end if
         
         if( ncycle .eq. 0 ) then
c
c           Set up grid structure, allocate storage ...
c
            call qreset(sizeq)
            nx(lf) = nxf
            ny(lf) = nyf
            h(lf)  = hf
            u(lf)    = imposs
            rhs(lf)  = imposs
            mres(lf) = qalloc(nx(lf)*ny(lf))
            du(lf)   = mres(lf)
            do 1000 l = lf - 1 , lc , -1
               nx(l) = (nx(l+1) - 1) / 2 + 1
               ny(l) = (ny(l+1) - 1) / 2 + 1
               h(l) = 2.0d0 * h(l+1)

               u(l)    = qalloc(nx(l)*ny(l))
               rhs(l)  = qalloc(nx(l)*ny(l))
               mres(l) = qalloc(nx(l)*ny(l))
               du(l)   = mres(l)
 1000       continue
            if( ltrace ) then
               call ivdump(nx,nl,'nx',0)
               call ivdump(ny,nl,'ny',0)
               call dvdump(h,nl,'h',0)
               call ivdump(u,nl,'u pointers',0)
               call ivdump(rhs,nl,'rhs pointers',0)
               call ivdump(mres,nl,'mres pointers',0)
               call ivdump(du,nl,'du pointers',0)

            end if
         end if
c
c        Cycle from fine to coarse ...
c
         do l = lf , lc , -1
            if( l .ne. lf  .or.  ncycle .eq. 0 ) then
               if(      l .eq. lc ) then
 100              continue
                     resnrm = relaxc2(q(u(l)),q(rhs(l)),
     *                                nx(l),ny(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg2v0: pre sweep',l,resnrm)
                     end if
                  if( resnrm .gt. epsic ) go to 100
               else if( l .eq. lf ) then
                  do sweep = 1 , preswf
                     resnrm = relaxc2(uf,rhsf,
     *                                nx(l),ny(l),h(l))
                     call cmres2(q(mres(l)),uf,rhsf,
     *                           nx(l),ny(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg2v0: pre sweep',l,resnrm)
                     end if
                  end do
               else 
                  do sweep = 1 , presw 
                     resnrm = relaxc2(q(u(l)),q(rhs(l)),
     *                                nx(l),ny(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg2v0: pre sweep',l,resnrm)
                     end if
                  end do
               end if
            end if
c
c           Set up coarse grid problem ...
c
            if(      l .eq. lf ) then
               call cmres2(q(mres(l)),uf,rhsf,nx(l),ny(l),h(l))
            else if( l .ne. lc ) then
               call cmres2(q(mres(l)),q(u(l)),q(rhs(l)),
     *                    nx(l),ny(l),h(l))
            end if
            if( l .ne. lc ) then
               call dmrshw(q(mres(l)),q(rhs(l-1)),
     *                     nx(l-1),ny(l-1))
               call dmls(q(u(l-1)),0.0d0,
     *                   nx(l-1),ny(l-1))
            end if

         end do 
c
c        Cycle from coarse to fine ...
c

         do l = lc + 1 , lf 
c
c           Interpolate correction ...
c
            call dmlint(q(u(l-1)),q(du(l)),nx(l-1),ny(l-1))
            if( l .eq. lf ) then
               call dmma(uf,q(du(l)),uf,nx(l),ny(l))
            else 
               call dmma(q(u(l)),q(du(l)),q(u(l)),nx(l),ny(l))
            end if

            if( l .eq. lf ) then
               do sweep = 1 , pstswf
                  resnrm = relaxc2(uf,rhsf,nx(l),ny(l),h(l))
                  if( ltrace ) then
                     call trace_res('mg2v0: post sweep',l,resnrm)
                  end if
               end do
            else 
               do sweep = 1 , pstsw 
                  resnrm = relaxc2(q(u(l)),q(rhs(l)),
     *                             nx(l),ny(l),h(l))
                  if( ltrace ) then
                     call trace_res('mg2v0: post sweep',l,resnrm)
                  end if
               end do
            end if

         end do

         ncycle = ncycle + 1
         mg2v0 = resnrm

         return

      end


c-----------------------------------------------------------------------
c
c     Basic 3d v-cycle multi grid solver for 
c
c     Laplacian[u(x,y,z)] = rhs(x,y,z)
c
c     on rectangular domain. 
c
c     This version allows nxf, nyf, nzy which are of form
c     2^p * <small number> + 1
c
c     where <small number> .le. maxn0 defined in the parameter 
c     statement below
c
c-----------------------------------------------------------------------

      double precision function mg3v0
     *   (uf,rhsf,nxf,nyf,nzf, hf, q, sizeq,
     *    preswf,pstswf,presw,pstsw,ncycle)

         implicit      none

         real*8        d3nrm2,    relaxc3
         integer       qalloc,    lilog2
         logical       lp2np1

         integer       nxf,       nyf,     nzf,      sizeq
         real*8        uf(nxf,nyf,nzf),    rhsf(nxf,nyf,nzf),
     *                 q(sizeq)
         real*8        hf
         integer       preswf,    pstswf,  presw,    pstsw,
     *                 ncycle
c
c        Multi-grid parameters, data structures and pointers 
c
         real*8        epsic
         parameter   ( epsic = 1.0d-10 )

         integer       lmin,      lmax
         parameter   ( lmin = 1,  lmax = 20 )

         integer       nx(lmin:lmax),      ny(lmin:lmax),
     *                 nz(lmin:lmax)
         real*8        h(lmin:lmax)

         integer       u(lmin:lmax),       rhs(lmin:lmax),
     *                 mres(lmin:lmax),    du(lmin:lmax)

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

         integer          p,      pceil
         common  / comp / p,      pceil
c
c        For parameter (nxf, nyf, nzf) checking 
c
         integer       checkfacts

         logical       badargs

         integer       nfact_x,   nx0
         integer       nfact_y,   ny0
         integer       nfact_z,   nz0

         integer       maxn0
         parameter   ( maxn0 = 17 )
c
c        Locals ...
c
         integer       lxf,       lyf,       lzf,
     *                 l,         lf,        lc,       nl,
     *                 sweep

         real*8        resnrm

         logical       ltrace 
         parameter   ( ltrace = .false. )

         save

         mg3v0 = -1.0d0
c
c        Parameter check ...
c

         if( ltrace ) then
            write(0,*) '>>> mg3v0: Parameter dump: '
            write(0,*) nxf, nyf, nzf, hf
            write(0,*) preswf, pstswf, presw, pstsw, ncycle
         end if

         badargs = .false.
         if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then
            write(*,*) 'mg3v0: Bad nxf value: ', nxf
            badargs = .true.
         end if
         if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then
            write(*,*) 'mg3v0: Bad nyf value: ', nyf
            badargs = .true.
         end if
         if( checkfacts(nzf - 1,maxn0,nz0,nfact_z) .lt. 0 ) then
            write(*,*) 'mg3v0: Bad nzf value: ', nzf
            badargs = .true.
         end if

         if( badargs) then
            return
         end if

         nl = min(nfact_x,nfact_y,nfact_z) + 1
         lc = 1
         lf = lc + nl - 1

         nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1
         ny0 = (nyf - 1) / (2 ** (nl - 1)) + 1
         nz0 = (nzf - 1) / (2 ** (nl - 1)) + 1

         if( ltrace ) then
            write(*,*) 'mg3v0: ', nl, ' levels.'
            write(*,*) 'mg3v0:  Coarse grid ', nx0, ' x ', ny0,
     *                 ' x ', nz0
         end if
         
         if( ncycle .eq. 0 ) then
c
c           Set up grid structure, allocate storage ...
c
            call qreset(sizeq)
            nx(lf) = nxf
            ny(lf) = nyf
            nz(lf) = nzf
            h(lf)  = hf
            u(lf)    = imposs
            rhs(lf)  = imposs
            mres(lf) = qalloc(nx(lf)*ny(lf)*nz(lf))
            du(lf)   = mres(lf)
            do 1000 l = lf - 1 , lc , -1
               nx(l) = (nx(l+1) - 1) / 2 + 1
               ny(l) = (ny(l+1) - 1) / 2 + 1
               nz(l) = (nz(l+1) - 1) / 2 + 1
               h(l) = 2.0d0 * h(l+1)

               u(l)    = qalloc(nx(l)*ny(l)*nz(l))
               rhs(l)  = qalloc(nx(l)*ny(l)*nz(l))
               mres(l) = qalloc(nx(l)*ny(l)*nz(l))
               du(l)   = mres(l)
 1000       continue
            if( ltrace ) then
               call ivdump(nx,nl,'nx',0)
               call ivdump(ny,nl,'ny',0)
               call ivdump(nz,nl,'nz',0)
               call dvdump(h,nl,'h',0)
               call ivdump(u,nl,'u pointers',0)
               call ivdump(rhs,nl,'rhs pointers',0)
               call ivdump(mres,nl,'mres pointers',0)
               call ivdump(du,nl,'du pointers',0)

            end if
         end if
c
c        Cycle from fine to coarse ...
c
         do l = lf , lc , -1
            if( l .ne. lf  .or.  ncycle .eq. 0 ) then
               if(      l .eq. lc ) then
 100              continue
                     resnrm = relaxc3(q(u(l)),q(rhs(l)),
     *                                nx(l),ny(l),nz(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg3v0: pre sweep',l,resnrm)
                     end if
                  if( resnrm .gt. epsic ) go to 100
               else if( l .eq. lf ) then
                  do sweep = 1 , preswf
                     resnrm = relaxc3(uf,rhsf,
     *                                nx(l),ny(l),nz(l),h(l))
                     call cmres3(q(mres(l)),uf,rhsf,
     *                           nx(l),ny(l),nz(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg3v0: pre sweep',l,resnrm)
                     end if
                  end do
               else 
                  do sweep = 1 , presw 
                     resnrm = relaxc3(q(u(l)),q(rhs(l)),
     *                                nx(l),ny(l),nz(l),h(l))
                     if( ltrace ) then
                        call trace_res('mg3v0: pre sweep',l,resnrm)
                     end if
                  end do
               end if
            end if
c
c           Set up coarse grid problem ...
c
            if(      l .eq. lf ) then
               call cmres3(q(mres(l)),uf,rhsf,nx(l),ny(l),nz(l),h(l))
            else if( l .ne. lc ) then
               call cmres3(q(mres(l)),q(u(l)),q(rhs(l)),
     *                    nx(l),ny(l),nz(l),h(l))
            end if
            if( l .ne. lc ) then
               call d3rshw(q(mres(l)),q(rhs(l-1)),
     *                     nx(l-1),ny(l-1),nz(l-1))
               call d3ls(q(u(l-1)),0.0d0,
     *                   nx(l-1),ny(l-1),nz(l-1))
            end if

         end do 
c
c        Cycle from coarse to fine ...
c

         do l = lc + 1 , lf 
c
c           Interpolate correction ...
c
            call d3lint(q(u(l-1)),q(du(l)),nx(l-1),ny(l-1),nz(l-1))
            if( l .eq. lf ) then
               call d33a(uf,q(du(l)),uf,nx(l),ny(l),nz(l))
            else 
               call d33a(q(u(l)),q(du(l)),q(u(l)),nx(l),ny(l),nz(l))
            end if

            if( l .eq. lf ) then
               do sweep = 1 , pstswf
                  resnrm = relaxc3(uf,rhsf,nx(l),ny(l),nz(l),h(l))
                  if( ltrace ) then
                     call trace_res('mg3v0: post sweep',l,resnrm)
                  end if
               end do
            else 
               do sweep = 1 , pstsw 
                  resnrm = relaxc3(q(u(l)),q(rhs(l)),
     *                             nx(l),ny(l),nz(l),h(l))
                  if( ltrace ) then
                     call trace_res('mg3v0: post sweep',l,resnrm)
                  end if
               end do
            end if

         end do

         ncycle = ncycle + 1
         mg3v0 = resnrm

         return

      end


c-----------------------------------------------------------------------
c
c     Resets pointer into main storage array.
c
c-----------------------------------------------------------------------

      subroutine qreset(ceil) 

         implicit         none

         integer          ceil

         integer          p,      pceil
         common  / comp / p,      pceil

         pceil = ceil
         p = 1

         return

      end

c-----------------------------------------------------------------------
c
c     Returns next available location in main storage array and updates 
c     pointer.
c
c-----------------------------------------------------------------------

      integer function qalloc(size) 

         integer          size

         integer          p,      pceil
         common  / comp / p,      pceil

         if( size .gt. 0 ) then
            if( p + size .le. pceil ) then
               qalloc = p
               p = p + size
            else 
               write(0,1000) int(size / 1.0d6 + 0.5d0),
     &                       int(pceil / 1.0d6 + 0.5d0)
1000           format(' >>> qalloc: Out of memory'/
     &                ' >>> qalloc: Requested block size:', i5,
     &                ' Mbytes'/
     &                ' >>> qalloc: Total available:     ', i5,
     &                ' Mbytes')
               stop
            end if
         else 
            write(0,*) '>>> qalloc: size ', size, ' <= 0.'
            qalloc = -1
         end if
         
         return

      end

c-----------------------------------------------------------------------
c
c     Returns number of words of memory used.
c
c-----------------------------------------------------------------------

      integer function qused()

         integer          p,      pceil
         common  / comp / p,      pceil

         qused = p

         return

      end
c-----------------------------------------------------------------------
c
c     Returns number of words of available memory.
c
c-----------------------------------------------------------------------

      integer function qleft()

         integer          p,      pceil
         common  / comp / p,      pceil

         qleft = pceil - p + 1

         return

      end

c-----------------------------------------------------------------------
c
c     Tracing routine
c
c-----------------------------------------------------------------------

      subroutine trace_res(string,l,resnrm)
   
         implicit      none

         character*(*) string 
         integer       l
         real*8        resnrm

         write(0,1000) string, l, resnrm
1000     format(a,' Level: ',i2,'  |res.|: ',1pe12.3)

         return

      end

c-----------------------------------------------------------------------
c
c     Does one-step 1d Point-GS relaxation and returns norm of
c     running ('instantaneous') residuals---checkerboard ordering.
c
c-----------------------------------------------------------------------
 
      double precision function relaxc1(u,rhs,nx,h)
 
         implicit       none
 
         integer        nx
         real*8         u(nx),             rhs(nx)
         real*8         h
 
         real*8         hm2,      jacel,   rresl
         integer        i,        pass  
 
         relaxc1 = 0.0d0
         hm2 = 1.0d0 / (h * h)
         jacel = -2.0d0 * hm2
 
         do pass = 1 , 2
            do i = pass + 1 , nx - 1 , 2
               rresl = hm2 * (u(i+1) + u(i-1) - 2.0d0*u(i)) -
     *                 rhs(i)
               u(i) = u(i) - rresl / jacel
               relaxc1 = relaxc1 + rresl * rresl
            end do
         end do
         relaxc1 = sqrt(relaxc1 / nx)
 
         return
 
      end
 
c-----------------------------------------------------------------------
c
c     Computes negative of 1d residuals ...
c
c-----------------------------------------------------------------------
 
      subroutine cmres1(mres,u,rhs,nx,h)
 
         implicit       none
 
         integer        nx
         real*8         mres(nx),         u(nx),      
     *                  rhs(nx)
 
         real*8         h
 
         integer        i
         real*8         hm2
 
         hm2 = 1.0d0 / (h * h)
 
         do i = 2 , nx - 1
            mres(i) = rhs(i) - 
     *                (hm2 * (u(i-1) - 2.0d0 * u(i) + u(i+1)))
         end do
         call dvlbs(mres,0.0d0,nx)
 
         return
 
      end


c-----------------------------------------------------------------------
c
c     Does one-step 2d Point-GS relaxation and returns norm of
c     running ('instantaneous') residuals---checkerboard ordering.
c
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
c     Computes negative of 2d residuals ...
c
c-----------------------------------------------------------------------
 
      subroutine cmres2(mres,u,rhs,nx,ny,h)
 
         implicit       none
 
         integer        nx,               ny
         real*8         mres(nx,ny),      u(nx,ny),   
     *                  rhs(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
               mres(i,j) = rhs(i,j) -
     *                       hm2 * (u(i+1,j) + u(i-1,j) +
     *                              u(i,j+1) + u(i,j-1) -
     *                              4.0d0*u(i,j))
            end do
         end do
         call dmlbs(mres,0.0d0,nx,ny)
 
         return
 
      end

c-----------------------------------------------------------------------
c
c     Does one-step 3d Point-GS relaxation and returns norm of
c     running ('instantaneous') residuals---checkerboard ordering.
c
c-----------------------------------------------------------------------
 
      double precision function relaxc3(u,rhs,nx,ny,nz,h)
 
         implicit       none
 
         integer        nx,              ny,            nz
         real*8         u(nx,ny,nz),     rhs(nx,ny,nz)
         real*8         h
 
         real*8         hm2,      mjelm1,   rresl
         integer        i,        j,        k,
     *                  isw,      ksw,      pass
 
         relaxc3 = 0.0d0
         hm2 = 1.0d0 / (h * h)
         mjelm1 =  h * h / 6.0d0
 
         ksw = 2
          do pass = 1 , 2
            isw = ksw
             do k = 2 , nz - 1
                do j = 2 , ny - 1
                   do i = isw + 1 , nx - 1, 2
                     rresl = hm2 * (u(i+1,j,k) + u(i-1,j,k) +
     *                              u(i,j+1,k) + u(i,j-1,k) +
     *                              u(i,j,k+1) + u(i,j,k-1) -
     *                              6.0d0 * u(i,j,k)) -
     *                       rhs(i,j,k)
                     u(i,j,k) = u(i,j,k) + mjelm1 *  rresl 
                     relaxc3 = relaxc3 + rresl * rresl
                  end do
                  isw = 3 - isw
                end do
            end do
            ksw = 3 - ksw
         end do
         relaxc3 = sqrt(relaxc3 / (nx * ny * nz))
 
         return
 
      end
 
c-----------------------------------------------------------------------
c
c     Computes negative of 3d residuals ...
c
c-----------------------------------------------------------------------
 
      subroutine cmres3(mres,u,rhs,nx,ny,nz,h)
 
         implicit       none
 
         integer        nx,               ny,           nz
         real*8         mres(nx,ny,nz),   u(nx,ny,nz),
     *                  rhs(nx,ny,nz)
 
         real*8         h
 
         integer        i,                j,            k
         real*8         hm2
 
         hm2 = 1.0d0 / (h * h)
 
          do k = 2 , nz - 1
             do j = 2 , ny - 1
                do i = 2 , nx - 1
                  mres(i,j,k) = rhs(i,j,k) - 
     *                          hm2 * (u(i+1,j,k) + u(i-1,j,k) +
     *                                 u(i,j+1,k) + u(i,j-1,k) +
     *                                 u(i,j,k+1) + u(i,j,k-1) -
     *                                   6.0d0*u(i,j,k))
               end do
            end do
         end do
         call d3lbs(mres,0.0d0,nx,ny,nz)
 
         return
 
      end

c-----------------------------------------------------------------------
c
c     Helper routine for argument checking ...
c 
c     This version doesn't demand that 
c 
c     num = prime * 2^n
c
c-----------------------------------------------------------------------

      integer function checkfacts(num,maxbase,base,npow2)

         implicit      none

         integer       ivprimedecomp
         
         integer       num,        maxbase,     base,       npow2

         integer       maxnfact
         parameter   ( maxnfact = 100 )
         integer       vfact(maxnfact)
         integer       nfact,      i 

         logical       ltrace
         parameter   ( ltrace = .false. )
         
         if( ivprimedecomp(num,vfact,nfact) .gt. 0 ) then
            npow2 = 0
            do i = 1 , nfact - 1
               if( vfact(i) .ne. 2 ) goto 100
               npow2 = npow2 + 1
            end do
100         continue
            base = 1
            do i = npow2 + 1 , nfact
               base = base * vfact(i)
            end do
            if( base .le. maxbase ) then
               checkfacts = npow2
            else 
               checkfacts = -1
            end if
         else 
            checkfacts = -1
         end if
         if( ltrace ) then
            write(*,*) 'checkfacts: Number being factored ', num
            call ivdump(vfact,nfact,'factors',6)
         end if

         return

      end

c-----------------------------------------------------------------------
c     Computes and dumps true residuals
c-----------------------------------------------------------------------
      subroutine cres1_trace(label,unit,res,u,rhs,nx,h)
         implicit       none

         character*(*)  label
         integer        nx,     unit
         real*8         res(*), u(*), rhs(*)
         real*8         h

         integer        shape(1)
         real*8         vh(1)

         shape(1) = nx
         vh(1)    = h

         call cres(res,u,rhs,1,shape,vh)
         write(unit,1000) label
1000     format('cres1_trace: ',a)
         call dvdump(u,nx,'u',unit)
         call dvdump(res,nx,'residual',unit)

         return
      end

c-----------------------------------------------------------------------
c     Computes true (non-negated) residual.
c-----------------------------------------------------------------------
      subroutine cres(res,u,rhs,rank,shape,vh)

         implicit     none

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

         call cmres(res,u,rhs,rank,shape,vh)

         if(      rank .eq. 1 ) then
            call dvsm(res,-1.0d0,res,shape(1))
         else if( rank .eq. 2 ) then
            call dvsm(res,-1.0d0,res,shape(1)*shape(2))
         else if( rank .eq. 3 ) then
            call dvsm(res,-1.0d0,res,shape(1)*shape(2)*shape(3))
         else
            write(0,*) 'cmres: rank = ', rank,' not implemented'
         end if

         return
      end
         
c-----------------------------------------------------------------------
c     Front end to cmres1, cmres2, cmres3 ...
c-----------------------------------------------------------------------
      subroutine cmres(mres,u,rhs,rank,shape,vh)

         implicit     none

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

         if(      rank .eq. 1 ) then
            call cmres1(mres,u,rhs,shape(1),                  vh(1))
         else if( rank .eq. 2 ) then
            call cmres2(mres,u,rhs,shape(1),shape(2),         vh(1))
         else if( rank .eq. 3 ) then
            call cmres3(mres,u,rhs,shape(1),shape(2),shape(3),vh(1))
         else 
            write(0,*) 'cmres: rank = ', rank,' not implemented'
         end if

         return

      end

c-----------------------------------------------------------------------
c     Eventually should be incorporated into dmatlib.f
c
c     Multiply matrix with scalar ...
c-----------------------------------------------------------------------
      subroutine dmsm(a1,s1,a2,d1,d2)

         implicit         none

         integer          d1,      d2

         real*8           a1(d1,d2),        a2(d1,d2)
         real*8           s1

         call dvsm(a1,s1,a2,d1*d2)

         return

      end