      program dfc
c ======================================================================
c
c
c        programmer     Jos van Kan
c        version 1.1    date   10-04-97
c          using module mg preset instead of common block cmgconst
c        version 1.0    date   01-01-97
c
c
c
c   copyright (c) 1997
c   all rights reserved. no part of this publication may be reproduced,
c   stored in a retrieval system ( e.g., in memory, disk, or core)
c   or be transmitted by any means, electronic, mechanical, photocopy,
c   recording, or otherwise, without written permission from the
c   author.
c
c   pilotcode 2d mgpackage
c
c ********************************************************************
c
c                         DESCRIPTION
c
c   Mg module test in defect correction mode.
c
c ********************************************************************
c
c                      KEYWORDS
c
c ********************************************************************
c
c                      MODULES USED
      use mgpreset
c     use testmod
c
c ********************************************************************
c
c                      INPUT / OUTPUT    PARAMETERS
c
      implicit none
c
c ********************************************************************
c
c                      COMMON BLOCKS
c
c ********************************************************************
c
c                      LOCAL PARAMETERS
      integer, parameter :: double = kind (1d0)
      integer matsize, vecsize, bufsize
      real (kind = double), dimension (:), allocatable :: mat, rhs, sol
      real (kind = double), dimension (:), allocatable :: res, cor, ap
      real (kind = double), dimension (:), allocatable :: rbuf
      integer ibuf (300)
      double precision phi
      integer nxf, nyf, bval, mindim
      integer iter, nsmooth, outrbuf, ilevel, hilevel
      double precision oldres, newres, startres, redfac
      integer maparr (9)
      data maparr /2, 3, 4, 5, 1, 6, 7, 8, 9/
c
c ********************************************************************
c
c                      SUBROUTINES CALLED
c
      double precision  maxvec
c   FILLMAT          fills matrix on finest level.
c   maxvec           calculates infinity norm of a vector
c   PRELUDE          preliminary activities for mg
c   PRINTMAT         prints matrix in stencil form
c
c     
c ********************************************************************
c
c                      I/O
c
c ********************************************************************
c
c                      ERROR MESSAGES
c
c ********************************************************************
c
c                      PSEUDO CODE
c
c ======================================================================
      write (*, *) 'Enter nxf, nyf, nsmooth, bval'
      read (*,*) nxf, nyf, nsmooth, bval
      hilevel = 0
      mindim = min (nxf, nyf)
      do while (mindim > 2)
        hilevel = hilevel + 1
        mindim = mindim/2 + 1
      end do
 
      write (*, *) nxf, nyf, nsmooth, bval
      print *, "hilevel: ", hilevel
c     Write (*, *) 'Enter nxf, nyf, hilevel, bval, tr, tp'
c     read (*,*) nxf, nyf, hilevel, bval, tr, tp
c     write (*, *) nxf, nyf, hilevel, bval, tr, tp
c      nsmooth = 2
      phi = 90d0
      vecsize = (nxf+1)*(nyf+1)
      matsize = 9*vecsize
      allocate (mat (matsize), rhs (vecsize), sol (vecsize))
      allocate (cor (vecsize), res (vecsize), ap (vecsize))
      call autosize (nxf, nyf, FULLCORS, hilevel, bufsize,
     .               0,1,0,1)
      write (*, '("bufferspace needed: ", I8)')
     .     bufsize
c      call pmap ( maparr )
      allocate (rbuf (bufsize)) 
      ibuf (slacki) = 0
      ibuf (slacki + 1) = 1
      ibuf (slacki + 2) = 0
      ibuf (slacki + 3) = 1

      call fillmat (mat, nxf, nyf, phi, bval, ibuf (slacki))
      call prelude ( mat, nxf, nyf, FULLCORS, hilevel, rbuf,
     .                bufsize,
     .               0, 1, 0, 1, ibuf)
      call fillvec (rhs, nxf, nyf, ibuf (slacki))
      call zeroize (sol, nxf, nyf, 
     .ibuf(slacki),ibuf(slacki+1),ibuf(slacki+2),ibuf(slacki+3))
      call copyvec (rhs, res, nxf, nyf,
     .ibuf(slacki),ibuf(slacki+1),ibuf(slacki+2),ibuf(slacki+3))
      oldres = maxvec (rhs, nxf, nyf, ibuf (slacki)) 
      startres = oldres
      iter = 0
100   continue
      call calcres ( mat, sol, rhs, res, nxf, nyf, ap,  
     .ibuf(slacki),ibuf(slacki+1),ibuf(slacki+2),ibuf(slacki+3))
      newres = maxvec (res, nxf, nyf,  
     .ibuf(slacki),ibuf(slacki+1),ibuf(slacki+2),ibuf(slacki+3))
      redfac = newres / oldres
      write (*, '("residual: ", d16.8, "  conv fac: ", f16.8)') newres,
     .   redfac
c      if ((newres/startres) .gt. 1d-8 .and. iter .lt. 20) then
      if ((newres/startres) .gt. 1d-10 .and. iter .lt. 25) then
c      if ((newres/startres) .gt. 1d-10) then
      call v1step ( mat, res, nsmooth, cor, rbuf, ibuf)
c      call v1stepn ( mat, res, nsmooth, cor, res, rbuf, ibuf)
c     call v1gs ( mat, res, nsmooth, cor, rbuf, ibuf)
c     call f1step ( mat, res, nsmooth, cor, rbuf, ibuf)
      call triadvec (sol, 1d0, cor, sol, nxf, nyf,  
     .ibuf(slacki),ibuf(slacki+1),ibuf(slacki+2),ibuf(slacki+3))
      oldres = newres
      iter = iter + 1
      go to 100
      endif 
      write (*, '("reduction: ", d16.8, " in ", i4, " mg iterations")')
     .      newres/startres, iter
      write (*, '("asymptot reduction factor per smoothing step :",
     .      f16.8)') redfac**(1d0/nsmooth)
      write (*, '("average  reduction factor per smoothing step :",
     .      f16.8)') (newres/startres)**(1d0/(iter*nsmooth))
      end
