program multigrid_1D
!
!  Iterative solution of -d^2u/dx^2=0
!
      implicit double precision(a-h,o-z)
      parameter (maxLev=10,maxNx=500)
      dimension u(maxLev,0:maxNx), rhs(maxLev,0:maxNx)
      dimension du(0:maxNx),resid(0:maxNx)
      dimension imax(maxLev)
      dimension nit(maxLev)
!
!     Read Data from screen
!     ---------------------
      open(1,file='mgDat.dat')
         read(1,*) nsegC        ! segments on the coarsest level      
         read(1,*) ncyc         ! multigrid cycles
         read(1,*) mgLev        ! multigrid Levels
         if (mgLev.gt.maxLev) stop 'Increase maxLev'
         do lev=1,mgLev
            read(1,*) nit(lev)
         enddo
         if (mgLev.eq.1) nit(1) = 1
         read (1,*) kmod1, coef1
         read (1,*) kmod2, coef2
         read (1,*) kmod3, coef3
      close(1)
!
!     Set-up grids
!     ------------
      nx = nsegC  ! since nodes numbering starts from 0
      do lev=1,mgLev
         if (nx.gt.maxNx) stop 'Increase maxNx'
         imax(lev) = nx
         write(*,*)'Level= ',lev,'0-',nx
         nx = 2*nx
      enddo
      write(*,*)' Nodes on finest grid = ',imax(mgLev) + 1
      pause
!
!     Initialization (only for the finest grid)
!     -----------------------------------------
      open(1,file='init.dat')
      pi  = 4.d0*datan(1.d0)
      lev = mgLev
      nx  = imax(lev)      
      do i=0,nx
         xi = dfloat(i)/dfloat(nx)
         a  =     xi*dfloat(kmod1)*coef1
         a  = a + xi*dfloat(kmod2)*coef2
         a  = a + xi*dfloat(kmod3)*coef3
         u  (lev,i) = dsin(a*pi) 
         rhs(lev,i) = 0.d0
         write(1,*) i,xi,u(lev,i)
      enddo
      close(1)
!
!     All other levels
!     ----------------
      do lev=1,mgLev-1
      do i=0,nx
         u  (lev,i) = 0.d0
         rhs(lev,i) = 0.d0
      enddo
      enddo
!
!
      open (1,file='conv.dat')
      cost  = 0.d0
!     ============
      do ic=1,ncyc
!     ============
         write(*,*)'----------------------------------------------'
!        Compute on all but the coarse levels
!        -----------------
         do lev=mgLev,1,-1
!        -----------------
            if (mgLev.gt.1) then
               if (lev.eq.1) exit
               if (lev.ne.mgLev) call zero(maxLev,maxNx,lev,imax,u)
            endif
!           Iterative solution 
            do it1=1,nit(lev)
               cost = cost + 1.d0/2.d0**(mgLev-lev)

               call solve(maxNx,maxLev,lev,imax,u,du,rhs,resid,restot)
               write(*,2) lev,dlog10(restot)
               if (lev.eq.mgLev) write(1,3)cost,dlog10(restot)
            enddo
!           Residual restriction to coarser grid
            if (lev.ne.1) call restriction(lev,lev-1,maxLev,maxNx,
     :                                   imax,u,rhs)
!        -----
         enddo
!        -----
! 
!        Compute on all but the fine levels
!        --------------
         do lev=1,mgLev
!        --------------
            if (lev.eq.mgLev) exit
            do it2=1,nit(lev)
               cost = cost + 1.d0/2.d0**(mgLev-lev)
               call solve(maxNx,maxLev,lev,imax,u,du,rhs,resid,restot)
               write(*,2)lev,dlog10(restot)
            enddo
!           Prolongation to finer grid
            if (lev.lt.mgLEv) 
     :          call prolongation(lev+1,lev,maxLev,maxNx,imax,u,rhs)
!        -----
         enddo
!        -----
!     =====
      enddo 
!     =====
      close(1)
  2   format('Level ',i3,'  Residual norm = ', f15.9)
  3   format(f10.4, 2(1x,f15.9))
!
      end 
!
!
!*********************************************************************** 
      subroutine zero(maxLev,maxNx,lev,imax,u)
!*********************************************************************** 
      implicit double precision(a-h,o-z)
      dimension u(maxLev,0:maxNx)
      dimension imax(maxLev)
!
      nx = imax(lev)
      do i=0,nx
         u(lev,i) = 0.d0
      enddo
!
      return
      end
!
!
!***********************************************************************
      subroutine solve(maxNx,maxLev,lev,imax,u,du,rhs,resid,
     :                 restot)
!***********************************************************************
      implicit double precision(a-h,o-z)
      dimension u(maxLev,0:maxNx), rhs(maxLev,0:maxNx)
      dimension du(0:maxNx),resid(0:maxNx)
      dimension imax(maxLev)
!      
      c1p = -1.d0       ! coef of i+1
      c1m = -1.d0       ! coef of i-1
      cc  =  2.d0       ! coef of i
      nx  = imax(lev)
      dx  = 1.d0/dfloat(nx)
!
!     zero du
      do i=0,nx
         du(i) = 0.d0
      enddo 
! 
      restot = 0.d0
      do i=1,nx-1
         au = -u(lev,i-1) + 2.d0*u(lev,i) - u(lev,i+1)
         resid(i) = rhs(lev,i) - au
         du(i) = (dx*dx*rhs(lev,i) - au - c1m*du(i-1) - c1p*du(i+1))/cc
         restot = restot + resid(i)*resid(i)
      enddo
      restot = dsqrt(restot)/dfloat(nx+1)
!
!     Update solution
      do i=1,nx-1
         u(lev,i)=u(lev,i) + du(i) 
      enddo
!      
      return
      end
!
!
!*********************************************************************** 
      subroutine restriction(levF,levC,maxLev,maxNx,imax,u,rhs)
!***********************************************************************       
! Restriction (fine-to-coarse)
!  injection: r^H(i) = r^h(2*i)      
!
      implicit double precision(a-h,o-z)
      dimension u(maxLev,0:maxNx), rhs(maxLev,0:maxNx)
      dimension imax(maxLev)
!
      nxC = imax(levC)   ! nodes on coarse grid
      nxF = imax(levF)   ! nodes on fine grid
      dxF = 1.d0/dfloat(nxF)
!
      do ic=1,nxC-1
         jf = 2*ic    ! fine level node
         auF = -u(levF,jf-1) + 2.d0*u(levF,jf) - u(levF,jf+1)
         residF = rhs(levF,jf) - auF/dxF/dxF
         rhs(levC,ic) = residF
      enddo
!
      return
      end
!
!
!*********************************************************************** 
      subroutine prolongation(levF,levC,maxLev,maxNx,imax,u,rhs)
!***********************************************************************       
! Prolongation (coarse-to-fine)
! 
      implicit double precision(a-h,o-z)
      dimension u(maxLev,0:maxNx), rhs(maxLev,0:maxNx)
      dimension imax(maxLev)
!
      nxC = imax(levC)
      do ic=1,nxC-1     ! ic: coarse grid node
         jf = 2*ic      ! jf: fine grid node
         eh   = u(levC,ic) 
         eh1  = 0.5d0*( u(levC,ic ) + u(levC,ic+1) )
         u(levF,jf)   = u(levF,jf)   + eh
         u(levF,jf+1) = u(levF,jf+1) + eh1

         if (ic.eq.1) then
            eh1m = 0.5d0*( u(levC,ic-1) + u(levC,ic) )
            u(levF,jf-1) = u(levF,jf-1) + eh1m
         endif

      enddo
!
      return
      end