      program bicg3
c ======================================================================
c
c
c        programmer     Jos van Kan
c        version 1.0    date   23-03-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  BICG implements Henk's Bi-CgStab algorithm with
c  v1step3 as preconditioner.
c
c ********************************************************************
c
c                      KEYWORDS
c
c ********************************************************************
c
c         MODULES USED
c
      use mgpreset
c
c **********************************************************************
c
c                      INPUT / OUTPUT    PARAMETERS
c
      implicit none
c
c ********************************************************************
c
c                      COMMON BLOCKS
c
c
c ********************************************************************
c
c                      LOCAL PARAMETERS
      integer, parameter :: double = kind(1d0)
      integer, parameter :: MB = 1024**2
c      real (dp) :: mat, rbuf
c      real (dp) :: y, t, v, p
      real (double), dimension (:), allocatable :: mat 
      real (double), dimension (:), allocatable :: rbuf
      real (double), dimension (:), allocatable :: rhs, sol, res
      real (double), dimension (:), allocatable :: y, t, v, p
      real (double) :: rho0, rho1, alfa, beta, omega
      integer :: ibuf (300)
      integer :: nxf, nyf, nzf, bv
      integer :: vecsize, matsize, bufsize, heapstat
      integer :: iter, nsmooth, hilevel, mindim
      real (double) :: 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
c                      SUBROUTINES CALLED
c
      double precision  maxvec3, inp3
c   FILLMAT3          fills matrix on finest level.
c   inp3              calculates inner product
c   maxvec3           calculates infinity norm of a vector
c   PRELUDE3          preliminary activities for mg
c   PRINTMAT3         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, bv'
      read (*,*) nxf, nyf, nzf, nsmooth, bv
      write (*, *) nxf, nyf, nzf, nsmooth, bv
c      write (*,*) "Enter e1, e2, e3"
c      read (*,*) e1, e2, e3
      e1 = 1d0
      e2 = 1d0
      e3 = 1d0
      mindim = min (nxf, nyf, nzf)
c      nsmooth = 2
      hilevel = 0
      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: ", 7*vecsize
      write (*,'(a, f8.3, a)') "total: ", (8d0*(matsize + bufsize + 
     .    7*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), 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
      allocate (t (vecsize), y(vecsize), v(vecsize), p(vecsize), stat =
     .          heapstat )
      if (heapstat .ne. 0) then
         write (*,*) "could not allocate bicg-stab vectors" 
         stop
      else
         write (*,*) "bicg-stab vectors allocated" 
      endif
      t = 0d0
      y = 0d0
      v = 0d0
      p = 0d0
cc
      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, bv, 
     . 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 zeroize3 (sol, 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
      rho0 = 1d0
      alfa = 1d0
      omega = 1d0
      call zeroize3 (v, nxf, nyf, nzf, 
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
      call zeroize3 (p, nxf, nyf, nzf, 
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
100   continue
      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
      if ((newres/startres) .gt. 1d-10 .and. iter .le. 10) then
        rho1 = inp3 (rhs, res, nxf, nyf, nzf, 
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        beta = rho1/rho0*alfa/omega
        call triadv3 (p,-omega,v,p,nxf,nyf,nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        call triadv3 (res,beta,p,p,nxf,nyf,nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        call v1step3 ( mat, p, nsmooth, y, rbuf, ibuf)
        call mmv3 ( mat,y,v,nxf,nyf,nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        alfa = rho1/inp3 ( rhs, v, nxf, nyf, nzf, 
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        call triadv3 (sol,alfa,y,sol,nxf,nyf,nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        call triadv3 (res,-alfa,v,res,nxf,nyf,nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        call v1step3 ( mat, res, nsmooth, y, rbuf, ibuf)
        call mmv3 ( mat, y, t, nxf, nyf, nzf, 
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        omega = inp3 (t, res, nxf, nyf, nzf, 
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))/
     .          inp3 (t, t, nxf, nyf, nzf, 
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        call triadv3 (sol,omega,y,sol,nxf,nyf,nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        call triadv3 (res,-omega,t,res,nxf,nyf, nzf,
     . ibuf (slacki), ibuf (slacki + 1), ibuf (slacki + 2),
     . ibuf (slacki+3), ibuf (slacki + 4), ibuf (slacki + 5))
        rho0 = rho1
        iter = iter + 1
        oldres = newres
        go to 100
      endif 
      write (*, '("reduction: ", d16.8, " in ", i4, 
     .      " bicg-stab iterations")')
     .      newres/startres, iter
      write (*, '("average  reduction factor per smoothing step :",
     .      f16.8)') (newres/startres)**(1d0/(iter*2*nsmooth))
      end
