!! ---------------------------------------------------------
!! bicg.f90
!! 
!! This inverts a non-symmetric matrix using the stabilized
!! biconjugate gradient method of  Van der Vorst.
!!
!! The matrix system is A*x = b, where A is implicitly defined
!! by the matrix-vector multiply defined in the subroutine
!! matvec.  Matvec must have the form :
!!
!!       SUBROUTINE matvec(q,mx,my,apk)
!!
!! where q in the vector and apk is the matrix vector product.
!!
!!
!! ------------------------------------------------------------

SUBROUTINE bicg(mx,my,b,x,matvec, iter, resid, info)
  IMPLICIT REAL (a-h,o-z)
  INCLUDE "parms.h"


  EXTERNAL matvec
  REAL norm_b
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: q,b,x,r_hat,Ax, &
       &r, p, x_kmh,r_kmh,Ap,Arkmh

  norm_b = dot(b,b,mx,my)
  IF (ABS(norm_b) == 0) THEN
     x(1:mx,1:my) = 0.0
     iter = 1
     resid = 0
     RETURN
  END IF

  CALL matvec(x,mx,my,Ax)
  r(1:mx,1:my) = b(1:mx,1:my) - Ax(1:mx,1:my)
  r_hat(1:mx,1:my) = r(1:mx,1:my)
  p(1:mx,1:my) = r(1:mx,1:my)

  k = 1
  BiCG_loop : DO
     CALL matvec(p,mx,my,Ap)
     a = dot(r,r_hat,mx,my)/dot(Ap,r_hat,mx,my)
     x_kmh(1:mx,1:my) = x(1:mx,1:my) + a*p(1:mx,1:my)
     r_kmh(1:mx,1:my) = r(1:mx,1:my) - a*Ap(1:mx,1:my)
     CALL matvec(r_kmh,mx,my,Arkmh)
     omega_k = dot(r_kmh,Arkmh,mx,my)/dot(Arkmh,Arkmh,mx,my)
     x(1:mx,1:my) = x_kmh(1:mx,1:my) + omega_k*r_kmh(1:mx,1:my)
     dr = dot(r,r_hat,mx,my)
     r(1:mx,1:my) = r_kmh(1:mx,1:my) - omega_k*Arkmh(1:mx,1:my)
     b_k = a/omega_k*dot(r,r_hat,mx,my)/dr
     p(1:mx,1:my) = r(1:mx,1:my) + b_k*(p(1:mx,1:my) - &
          & omega_k*Ap(1:mx,1:my))

     !! Is residual small enough?
     IF (SQRT(dot(r,r,mx,my)) < resid) EXIT

     k = k + 1
     IF (k > iter) EXIT
  END DO BiCG_loop
  IF (k == iter) info = 1
  iter = k
  resid = SQRT(dot(r,r,mx,my))
  RETURN
END SUBROUTINE bicg


!! Dot product
REAL FUNCTION dot(A,B,mx,my)
  IMPLICIT REAL(a-h,o-z)

  INCLUDE "parms.h"

  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: A,B

  sum = 0.0
  DO i = 1,mx
     DO j = 1,my
        sum = sum + A(i,j)*B(i,j)
     END DO
  END DO
  dot = sum
END FUNCTION dot
