      program defect correction
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 3d 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, parameter :: MB = 1024**2
      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, nzf, bval
      integer iter, nsmooth, ilevel, hilevel, mindim
      integer :: vecsize, matsize, bufsize, heapstat
      double precision oldres, newres, startres, redfac
      real (double) :: e1, e2, e3
      integer :: maparr (27)
      data maparr /2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 14, 15,
     .      16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27/

c ********************************************************************
c
c                      SUBROUTINES CALLED
c
      double precision  maxvec3
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, nzf, nsmooth, bval'
      read (*,*) nxf, nyf, nzf, nsmooth, bval
      write (*, *) nxf, nyf, nzf, nsmooth, bval
c      nsmooth = 2
c      phi = 90d0
      e1 = 1d0
      e2 = 1d0
      e3 = 1d0
      mindim = min (nxf, nyf, nzf)
      hilevel = 0
      nsmooth = 2
      do while (mindim .gt. 2)
        hilevel = hilevel + 1
        mindim = mindim/2+1
      end do
      write (*,*) "hilevel: ", hilevel
      call autosiz3 (nxf, nyf, nzf, FULLCORS, hilevel, bufsize,
     .               0,1,0,1,0,1)
      write (*,'(a, i9)') "buffer space needed: ", bufsize
      matsize = (nxf+1)*(nyf+1)*(nzf+1)*27
      vecsize = matsize/27

      write (*,'(a, i9)') "matrix space needed: ", matsize
      write (*,'(a, i9)') "vector space needed: ", 5*vecsize
      write (*,'(a, f8.3, a)') "total: ", (8d0*(matsize + bufsize +
     .    5*vecsize))/MB, " MB allocatable space"

      allocate (mat (matsize), stat = heapstat)
      if (heapstat .ne. 0) then
         write (*,*) "could not allocate matrix, status = " , heapstat
         stop
      else
         write (*,*) "matrix allocated"
      endif
      mat = 0d0

      allocate (rbuf (bufsize), stat = heapstat)
      if (heapstat .ne. 0) then
         write (*,*) "could not allocate buffer, status= " , heapstat
         stop
      else
         write (*,*) "buffer allocated"
      endif
      rbuf = 0d0

      allocate (rhs (vecsize), sol (vecsize), res(vecsize), 
     .           cor (vecsize), ap (vecsize), stat = heapstat )
      if (heapstat .ne. 0) then
         write (*,*) "could not allocate problem vectors, status = " ,
     .               heapstat
         stop
      else
         write (*,*) "problem vectors allocated"
      endif

      rhs = 0d0
      sol = 0d0
      res = 0d0
      cor = 0d0
      ap = 0d0

      call pmap3 ( maparr )
      ibuf (slacki) = 0
      ibuf (slacki + 1) = 1
      ibuf (slacki + 2) = 0
      ibuf (slacki + 3) = 1
      ibuf (slacki + 4) = 0
      ibuf (slacki + 5) = 1

      call fillmat3 (mat, nxf, nyf, nzf, e1, e2, e3, bval,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
      call prelude3 ( mat, nxf, nyf, nzf, FULLCORS, hilevel, rbuf,
     .                bufsize,
     .               0, 1, 0, 1, 0, 1, ibuf)
      call fillvec3 (rhs, nxf, nyf, nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
      call copyvec3 (rhs, res, nxf, nyf, nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))

      oldres = maxvec3 (rhs, nxf, nyf, nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))

      startres = oldres
      iter = 0
100   continue
      call calcres3 ( mat, sol, rhs, res, nxf, nyf, nzf, ap,  
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
      newres = maxvec3 (res, nxf, nyf, nzf,  
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
      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. 15) then
c      if ((newres/startres) .gt. 1d-10) then
      call v1step3 ( 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 triadv3 (sol, 1d0, cor, sol, nxf, nyf, nzf,  
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
      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
