      program bicg
c ======================================================================
c
c        programmer     Jos van Kan
c        version 1.2    date   07-03-1997
c        version 1.0    date   01-01-1997
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  BICG implements Henk's Bi-CgStab algorithm with
c  v1stepn as preconditioner.
c **********************************************************************
c
c                       KEYWORDS
c
c **********************************************************************
c
c                       MODULES USED
c
      use mgpreset
c      use testmod
c **********************************************************************
c
c                       INPUT / OUTPUT PARAMETERS
c
      implicit none

c ********************************************************************
c                      COMMON BLOCKS
c
c **********************************************************************
c
c                       LOCAL PARAMETERS
c
      integer, parameter :: dp = kind (1d0)
      real (kind = dp), dimension (:), allocatable :: mat, rhs, sol
      real (kind = dp), dimension (:), allocatable :: res, y, t
      real (kind = dp), dimension (:), allocatable :: p, v
      real (kind = dp) :: rho0, rho1, alfa, beta, omega

      real (kind = dp), dimension (:), allocatable :: rbuf
      integer, dimension (300) :: ibuf 
      real (kind = dp) :: phi
      integer :: nxf, nyf, matsize, vecsize, bufsize, status
      integer :: iter, nsmooth, hilevel, bval, mindim

      real (kind = dp) :: oldres, newres, startres, redfac
      integer :: maparr (9)
      data maparr /2, 3, 4, 5, 1, 6, 7, 8, 9/

      real (kind = dp) ::  maxvec, inp

c   FILLMAT          fills matrix on finest level.
c   inp              calculates inner product
c   maxvec           calculates infinity norm of a vector
c   PRELUDE          preliminary activities for mg
c   PRINTMAT         prints matrix in stencil form
c **********************************************************************
c
c                       I/O
c
c **********************************************************************
c
c                       ERROR MESSAGES
c
c **********************************************************************
c
c                       PSEUDO CODE
c
c **********************************************************************
c
c                       COMMON BLOCKS
c
c **********************************************************************
c
c                       SUBROUTINES CALLED
c
c **********************************************************************
c
c                       DATA STATEMENTS
c
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

      call autosize (nxf, nyf, FULLCORS, hilevel, bufsize,
     +               0,1,0,1)
      write (*, *) 'Bufferspace needed: ', bufsize
      vecsize = (nxf + 1)*(nyf+1)
      matsize = 9 * vecsize
      allocate (mat (matsize), rhs (vecsize), sol (vecsize), 
     +           res(vecsize), p (vecsize), t(vecsize), 
     +           y (vecsize), v (vecsize), stat = status)
      write (*,*) "Space needed by matrix & vectors", 16*vecsize
      if (status .ne.0) then
      write (*,*) 'error in allocation of matrix/vectors'
      stop
      endif
      allocate (rbuf (bufsize), stat = status)
      if (status .ne.0) then
      write (*,*) 'error in allocation of rbuf'
      stop
      endif
      call pmap ( maparr )
      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)

c         write (*,*) "preparation time :", preptime, " seconds"

      call fillvec (rhs, nxf, nyf, ibuf (slacki))

c     --- tottime = 0e0
c         do itiming = 1, 10

      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), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))

c     --- elapsed_start = SECNDS (0e0)
c         ittime = 0e0

      startres = oldres
      iter = 0
      rho0 = 1d0
      alfa = 1d0
      omega = 1d0
      call zeroize (v, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
      call zeroize (p, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
100   continue
      newres = maxvec (res, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
      redfac = newres / oldres

c     --- if ((newres/startres).lt.1d-8) then
      write (*, '("residual: ", d16.8, "  conv fac: ", f16.8)') newres,
     +   redfac

c      if ((newres/startres).gt.1d-10 .and. iter.le.15) then
      if ((newres/startres).gt.1d-8 .and. iter.le.15) then
        rho1 = inp (rhs, res, nxf, nyf, 
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        beta = rho1/rho0*alfa/omega
        call triadvec (p, -omega, v, p, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        call triadvec (res, beta, p, p, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
c       call v1stepn ( mat, p, nsmooth, y, v, rbuf, ibuf)
       call v1step ( mat, p, nsmooth, y, rbuf, ibuf)
c        call v1gs ( mat, p, nsmooth, y, rbuf, ibuf)
c       call triadvec (p, -1d0, v, v, nxf, nyf,
c     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
c     + ibuf (slacki + 3))
c --- v now contains Ay
        call multmatv ( mat, y, v, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        alfa = rho1/inp ( rhs, v, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        call triadvec (sol, alfa, y, sol, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        call triadvec (res, -alfa, v, res, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
c  --- v1stepn gives the new residual on exit in t, which saves a 
c      matrix vector multiply
c       call v1stepn ( mat, res, nsmooth, y, t, rbuf, ibuf)
        call v1step ( mat, res, nsmooth, y, rbuf, ibuf)
c        call v1gs ( mat, res, nsmooth, y, rbuf, ibuf)
c       call triadvec (res, -1d0, t, t, nxf, nyf,
        call multmatv ( mat, y, t, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
c--- t now contains Ay
        omega = inp (t, res, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))/
     +          inp (t, t, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        call triadvec (sol, omega, y, sol, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        call triadvec (res, -omega, t, res, nxf, nyf,
     + ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     + ibuf (slacki + 3))
        rho0 = rho1
        iter = iter + 1
        oldres = newres
        go to 100
      end if


      write (*, '("reduction: ", d16.8, " in ", i4,
     +      " bicg-stab iterations")')
     +      newres/startres, iter
      write (*, '("asymptot reduction factor per smoothing step :",
     +      f16.8)') redfac**(1d0/(2*nsmooth))
      write (*, '("average  reduction factor per smoothing step :",
     +      f16.8)') (newres/startres)**(1d0/(iter*2*nsmooth))


      end
