*
*
************************************************************************
*

      subroutine bicgst (ilower,iupper,jlower,jupper,
     &             cc,cn,cs,ce,cw,u,rhs,
     &             eps,rmax,ier)

*
************************************************************************
*
*     This routine is a wrapper for the BI-CGStab solver.  The input
*  arrays are the finite difference coefficients on the 2D grid with
*  the boundary conditions absorbed into them.
*
****  Parameters
*
*  INTEGERS:
*     im,jm :
*          These define sizes of the coefficient arrays passed to this
*          routine.  These values are provided via the global header
*          file param.h.
*
*     ilower,iupper,
*     jlower,jupper :
*          These are the indices of the computational grid that
*          correspond to upper and lower index limits for points on
*          which computations are actually done.
*
*     ier :
*          This is an error flag with the following possible return
*          values:
*              ier =  0 => solver converged without error
*              ier = -1 => solver did not converge
*
*  REALS:
*     cc(im,jm),cn(im,jm),
*     cs(im,jm),ce(im,jm),cw(im,jm):
*          These are the finite difference coefficient arrays for a
*          five-point stencil on the two dimensional grid as follows:
*
*
*                      cn       
*            ^
*            |  cw     cc    ce
* increasing |
*  j values  |         cs       
*   (theta)  |
*             --------> increasing i values (eta)
*
*  Ie. the coresspondence is : n  : i,j+1
*                              ne : i+1,j+1
*                              c  : i,j
*                              e  : i+1,j
*                              s  : i,j-1
*
*     u(im,jm) :
*          Input: this contains the initial guess to the solution of the
*                 equation
*          Output: This contains the final approximation to the solution
*                 determined by the BI-CGStab solver.
*     
*     rhs(im,jm) :
*          This array contains the values of the right hand side of the
*          equation at every point on the rwo dimensional grid.
*
*     eps :
*          eps > 0.  The maximum norm of the residual is calculated 
*                    each BI-CGStab iteration. The algorithm is
*                    terminated when this maximum becomes less than eps 
*                    or when the maximum number of iterations (see ncyc)
*                    is exceeded.  It is up to the user to provide a
*                    meaningfull tolerance criteria for the particular
*                    problem being solved.
*     rmax:
*          This is the final value of the residual calculated.
*
************************************************************************
*
*
      implicit none

      include 'param.h'

      integer ilower,iupper,jlower,jupper,ier
      real cc(im,jm),cn(im,jm),cs(im,jm),ce(im,jm),
     &     cw(im,jm),u(im,jm),rhs(im,jm)
      real eps,rmax

*
************************************************************************
*
*  Variable definitions:
*
*  INTEGERS:
*
*     ncyc :
*          The maximum number of BI-CGStab cycles to be used. If
*          the maximum norm of the residual is not less than eps at the
*          end of ncyc cycles, the algorithm is terminated.
*          
*     idim,jdim :
*          The number of grid points in the x- and y-directions
*          respectively (including two ficticious points).  These
*          values are provided via the gloabal header file param.h.
*
*     iscale :
*          Flag to indicate whether problem can be diagonally scaled to
*          speed convergence of the multigrid solver.
*
*     i,j:
*          Loop/counting variables.
*
*  REALS:
*
*     a(idim,jdim,5) :
*         Contains all the finite difference coefficient values to be
*         passed to bicgstab as follows.
*               a(idim,jdim,1) = cc(idim+ilower,jdim+jlower)
*               a(idim,jdim,2) = cn(idim+ilower,jdim+jlower)
*               a(idim,jdim,3) = cs(idim+ilower,jdim+jlower)
*               a(idim,jdim,4) = ce(idim+ilower,jdim+jlower)
*               a(idim,jdim,5) = cw(idim+ilower,jdim+jlower)
*
*     f(idim,jdim) :
*          f contains the right hand side vector of the matrix
*          equation to be solved by bicgstab.
*
*     q(idim,jdim) :
*          q contains the initial guess on the fine grid.
*
************************************************************************
*
*

      integer ncyc
      parameter (ncyc=idim*jdim)

      real a(idim,jdim,5),
     &     q(idim,jdim),f(idim,jdim)

      integer iscale,i,j

*
************************************************************************
*

      ier=0

*
*  Set up coefficients into single array with correct indexing
*
      do 100 j=1,jdim
      do 100 i=1,idim
      a(i,j,1)=cc(i+ilower-1,j+jlower-1)
      a(i,j,2)=cn(i+ilower-1,j+jlower-1)
      a(i,j,3)=cs(i+ilower-1,j+jlower-1)
      a(i,j,4)=ce(i+ilower-1,j+jlower-1)
      a(i,j,5)=cw(i+ilower-1,j+jlower-1)
      q(i,j)=u(i+ilower-1,j+jlower-1)
      f(i,j)=rhs(i+ilower-1,j+jlower-1)
 100  continue



*
*  Determine whether we can diagonally scale the problem to speed
*  convergence. Can only be done if there are no zeros on the main
*  diagonal (ie. central difference coefficient).
*
      iscale=1
      do 200 j=1,jdim
      do 200 i=1,idim
      if (a(i,j,1) .eq. 0.) then
        iscale=0
      endif
 200  continue

*
*  Do the diagonal scaling if we can.
*
      if (iscale.eq.1) then

      do 210 j=1,jdim
      do 210 i=1,idim
        f(i,j)=f(i,j)/a(i,j,1)
        a(i,j,5)=a(i,j,5)/a(i,j,1)
        a(i,j,4)=a(i,j,4)/a(i,j,1)
        a(i,j,3)=a(i,j,3)/a(i,j,1)
        a(i,j,2)=a(i,j,2)/a(i,j,1)
        a(i,j,1)=1.
 210  continue

      endif

* 
*   Now call the bicgstab routine
*
      
      call bicgstab (a,q,f,eps,ncyc,rmax,ier)

      if (rmax.gt.eps) then
        ier = -1
        print*,'Did not converge'
        print*,' maximum residual    = ',rmax
        print*,' tolerance           = ',eps
      endif

*
*  Put solution into user's array.
*
      do 500 j=1,jdim
      do 500 i=1,idim
      u(i+ilower-1,j+jlower-1)=q(i,j)
 500  continue

      return
      end
*
*
************************************************************************
*

      subroutine bicgstab (a,x,r,tol,ncyc,rnorm,ier)

*
************************************************************************
*
*     This routine performs the BI-CGStab algorithm to solve the system
* of linear equations.
*
************************************************************************
*
****  Variable Definitions:
*
*  INTEGERS:
*
*     idim,jdim:
*          The number of points in the x- and y-directions respectively
*          on which computations are performed.
*
*     ncyc :
*          The maximum number of BI-CGStab cycles to be used. If
*          the maximum norm of the residual is not less than eps at the
*          end of ncyc cycles, the algorithm is terminated.
*
*     ier :
*          This is an error flag with the following possible return
*          values:
*              ier =  0 => solver converged without error
*              ier = -1 => solver did not converge
*
*  REALS:
*
*     a(idim,jdim,5):
*          Contains all the finite difference coefficients in a single
*          array for passing to the matvec routine.
*
*     x(idim,jdim):
*          Contains initial guess on input and final approximation to
*          the solution on return only in the region to be computed.
*
*     r(idim,jdim):
*          Contains the right-hand-side (source terms) of equations
*          only in the region to be computed.
*
*     tol :
*          tol > 0.  The maximum norm of the residual is calculated
*                    each BI-CGStab iteration. The algorithm is
*                    terminated when this maximum becomes less than tol
*                    or when the maximum number of iterations (see ncyc)
*                    is exceeded.  It is up to the user to provide a
*                    meaningfull tolerance criteria for the particular
*                    problem being solved.
*
*     rnorm:
*          This is the value of the maximum residual calculated.
*
***********************************************************************
*

      implicit none

      include 'param.h'

      integer ncyc,ier
      real a(idim,jdim,5),x(idim,jdim),r(idim,jdim),
     .     tol,rnorm


*
***********************************************************************
*
*
*  BiCGstab Algorithm Variables:
*
*  INTEGERS:
*
*     i,j,k:
*          Looping/counting variables.
*
*  REALS:
*
*     omega, chi, delta, deltap, phi:
*          Direction vector coefficients.
*
*     p(idim,jdim), Ap(idim,jdim), w(idim,jdim), As(idim,jdim):
*          Direction vectors and matrix vector product temporaries.
*
***********************************************************************
*

      integer i,j,k
      real omega, chi,chi1,chi2, delta, deltap, phi
      real p(idim,jdim), Ap(idim,jdim), w(idim,jdim), As(idim,jdim)

*
***********************************************************************
*

      k = 0

  10  call usermv(a,x,Ap)

      do 1000 j=1,jdim
      do 1000 i=1,idim
      r(i,j)=r(i,j)-Ap(i,j)
 1000 p(i,j)=r(i,j)

      delta = 0.
      do 100 j=1,jdim
      do 100 i=1,idim
 100  delta=delta+r(i,j)

      if (delta.eq.0.) then
        ier=-1
        return
      endif

      call usermv(a,p,Ap)

      phi = 0.
      do 110 j=1,jdim
      do 110 i=1,idim
 110  phi=phi+Ap(i,j)
      phi=phi/delta

      if (phi.eq.0.) then
        ier=-1
        return
      endif

      rnorm=0.
      do 120 j=1,jdim
      do 120 i=1,idim
 120   rnorm=rnorm+r(i,j)*r(i,j)
      rnorm=sqrt(rnorm)

c     Test if initial guess is great (residual less than tolerance)
      if (rnorm .lt. tol)  return

 1    continue
      k = k + 1

      omega = 1./phi

      do 1100 j=1,jdim
      do 1100 i=1,idim
 1100 w(i,j) = r(i,j) - omega*Ap(i,j)

      call usermv(a,w,As)

      chi1=0.
      chi2=0.
      do 130 j=1,jdim
      do 130 i=1,idim
      chi1=chi1+As(i,j)*w(i,j)
 130  chi2=chi2+As(i,j)*As(i,j)

      chi = chi1/chi2

      do 1200 j=1,jdim
      do 1200 i=1,idim
      r(i,j) = w(i,j) - chi*As(i,j)
 1200 x(i,j) = x(i,j) + omega*p(i,j) +chi*w(i,j)

      deltap = delta
      delta = 0.
      do 140 j=1,jdim
      do 140 i=1,idim
 140  delta=delta+r(i,j)

      if (delta.eq.0.) then
        goto 10
      endif

      do 1300 j=1,jdim
      do 1300 i=1,idim
 1300 p(i,j) = r(i,j) + (p(i,j) - chi*Ap(i,j))
     .                   *(omega*delta/(deltap*chi))

      call usermv(a,p,Ap)

      phi = 0.
      do 150 j=1,jdim
      do 150 i=1,idim
 150  phi=phi+Ap(i,j)
      phi=phi/delta

      if (phi.eq.0.) then
        goto 10
      endif

      if (k.gt.ncyc) then
        print*,' BI-CGStab solver reached maximum nuber of iterations.'
        ier=-1
        return
      endif

      rnorm=0.
      do 160 j=1,jdim
      do 160 i=1,idim
 160   rnorm=rnorm+r(i,j)*r(i,j)
      rnorm=sqrt(rnorm)


      if (rnorm .gt. tol) goto 1

      return
      end
*
*
***********************************************************************
*

      subroutine usermv(a,x,y)

*
***********************************************************************
*
*     This routine performs the matrix-vector multiply (A*x=y)
*  necessary for a 5-point finite difference stencil in 2D.  If the
*  parameter nproc>1, and the code is compiled for parallelization, the
*  matrix-vector multiply will be parallelized ideally for nproc
*  processors.
*
***********************************************************************
*
*  Variable Definitions
*
*  INTEGERS:
*
*     idim,jdim:
*          The number of points in the x- and y-directions respectively
*          on which computations are performed.
*
*     nproc:
*          When compiling for parallel executiuon, this is the number
*          of processors the code is intended to be run on.
*
*     i,j,ii:
*          Indexing/looping variables.
*
*  REALS:
*
*     a(idim,jdim,5):
*          The finite difference coefficient array.
*
*     x(idim,jdim),y(idim,jdim):
*          The source and result vectors of the matrix product A*x=y.
*
***********************************************************************
*

      implicit none

      include 'param.h'

      integer i,j,ii
      real a(idim,jdim,5),x(idim,jdim),y(idim,jdim)

      integer nproc
      parameter (nproc=1)

*
***********************************************************************
*
*

      i=1
      y(i,1) = a(i,1,1)*x(i,1)
     .        +a(i,1,2)*x(i+idim,1)
     .        +a(i,1,4)*x(i+1,1)

      do 150 i=2,idim
 150  y(i,1) = a(i,1,1)*x(i,1)
     .        +a(i,1,2)*x(i+idim,1)
     .        +a(i,1,4)*x(i+1,1)
     .        +a(i,1,5)*x(i-1,1)

      if (nproc.eq.1) then

      do 200 i=idim+1,idim*jdim-idim

 200  y(i,1) = a(i,1,1)*x(i,1)
     .        +a(i,1,2)*x(i+idim,1)
     .        +a(i,1,3)*x(i-idim,1)
     .        +a(i,1,4)*x(i+1,1)
     .        +a(i,1,5)*x(i-1,1)

      else

      ii = (idim*jdim-2*idim)/nproc

      do 210 j=1,nproc
      do 210 i=idim+(j-1)*ii+1,idim+j*ii
 210  y(i,1) = a(i,1,1)*x(i,1)
     .        +a(i,1,2)*x(i+idim,1)
     .        +a(i,1,3)*x(i-idim,1)
     .        +a(i,1,4)*x(i+1,1)
     .        +a(i,1,5)*x(i-1,1)

      do 211 i=idim+nproc*ii+1,idim*jdim-idim
 211  y(i,1) = a(i,1,1)*x(i,1)
     .        +a(i,1,2)*x(i+idim,1)
     .        +a(i,1,3)*x(i-idim,1)
     .        +a(i,1,4)*x(i+1,1)
     .        +a(i,1,5)*x(i-1,1)

      endif

      do 250 i=idim*jdim-idim+1,idim*jdim-1

 250  y(i,1) = a(i,1,1)*x(i,1)
     .        +a(i,1,3)*x(i-idim,1)
     .        +a(i,1,4)*x(i+1,1)
     .        +a(i,1,5)*x(i-1,1)

      i=idim*jdim
      y(i,1) = a(i,1,1)*x(i,1)
     .        +a(i,1,3)*x(i-idim,1)
     .        +a(i,1,5)*x(i-1,1)

      return
      end
*
*
*
***********************************************************************
*
