!=======================================================================
      subroutine cgnr9(nx,ny,nbw,itmax,itCGNR,level,ierr,tol,
     &              A,AT,ILU,ILUT,x,b,wksp)
!=======================================================================
! Seongjai Kim,  skim@ms.uky.edu, Tue Mar 12 08:32:17 EST 2002
!=======================================================================
      implicit none
      integer nx,ny,nbw,itmax,itCGNR,level,ierr
      real*8  tol
      complex*16 A(-nbw:nbw,0:nx,0:ny),AT(-nbw:nbw,0:nx,0:ny)
      complex*16 ILU(-nbw:nbw,0:nx,0:ny),ILUT(-nbw:nbw,0:nx,0:ny)
      complex*16 x(0:nx,0:ny),b(0:nx,0:ny),wksp(0:nx,0:ny,4)

      integer m,ix,iy
      integer idp,idq,idr,idy,idz,neqn
      integer level0,level1
      real*8  rho0,rho1,r0,rm,alpha,beta
      complex*16 dotProd,ctmp

!-----------------------

      idp=1  ! wksp(,,idp) saves vector p
      idq=2  ! wksp(,,idq) saves vector Ap
      idr=3  ! wksp(,,idr) saves vector r
      idy=4  ! wksp(,,idy) saves vector A^*(Ap)
      idz=5  ! wksp(,,idz) saves vector z

      neqn=(nx+1)*(ny+1)
      level0=0
      level1=1

!-------------------------
! Print User Information
!-------------------------

      if(level.ge.1) then
          print'("CGNR9: nx=",i4," ny=",i4)',nx,ny
          print*,"itmax=",itmax," tol=",tol
      endif

      if(nbw.ne.4)then
          print*,"ERROR: cgnr9.f: Wrong nbw=",nbw
          ierr=1
          return
      endif


!---------------------------------------
! Get "AT", "ILU", and "ILUT" (:=ILU^*)
!---------------------------------------

      call ilu9(nx,ny,nbw,level,ierr,A,ILU)
      call transConjg(nx,ny,nbw,level,ierr,A,AT)
      call transConjg(nx,ny,nbw,level,ierr,ILU,ILUT)


!--------------------------------
! Initial setup
!--------------------------------

      call mtx9vec(nx,ny,A,x,wksp(0,0,idq),level,ierr)

      do iy=0,ny
      do ix=0,nx
          wksp(ix,iy,idy)=b(ix,iy)-wksp(ix,iy,idq)
      enddo
      enddo

      call mtx9vec(nx,ny,AT,wksp(0,0,idy),wksp(0,0,idr),level,ierr)

      ctmp=dotProd(neqn,level0,ierr,wksp(0,0,idr),wksp(0,0,idr))
      r0=dsqrt(dreal(ctmp))

      do iy=0,ny
      do ix=0,nx
          wksp(ix,iy,idz)=wksp(ix,iy,idr)
      enddo
      enddo

      call ilu9Subst(2,nx,ny,nbw,level0,ierr,ILUT,wksp(0,0,idz))
      call ilu9Subst(1,nx,ny,nbw,level0,ierr,ILU,wksp(0,0,idz))

      do iy=0,ny
      do ix=0,nx
          wksp(ix,iy,idp)=wksp(ix,iy,idz)
      enddo
      enddo

      ctmp=dotProd(neqn,level0,ierr,wksp(0,0,idz),wksp(0,0,idr))
      rho0=dreal(ctmp)

      if(r0.le.tol)then
         if(level.ge.1) print*,"X0 is accurate enough: r0=",r0
         return
      endif


!--------------------------------
! Begin The PCGNR Iteration
!--------------------------------
      do m=1,itmax
          
          call mtx9vec(nx,ny,A,wksp(0,0,idp),wksp(0,0,idq),level,ierr)
          ctmp=dotProd(neqn,level0,ierr,wksp(0,0,idq),wksp(0,0,idq))
          alpha=rho0/dreal(ctmp)

          do iy=0,ny
          do ix=0,nx
              x(ix,iy)=x(ix,iy)+alpha*wksp(ix,iy,idp)
          enddo
          enddo

          if(mod(m,100).eq.0)then
              call mtx9vec(nx,ny,A,x,wksp(0,0,idq),level,ierr)
              do iy=0,ny
              do ix=0,nx
                  wksp(ix,iy,idy)=b(ix,iy)-wksp(ix,iy,idq)
              enddo
              enddo
              call mtx9vec(nx,ny,AT,wksp(0,0,idy),wksp(0,0,idr),0,ierr)
          else
              call mtx9vec(nx,ny,AT,wksp(0,0,idq),wksp(0,0,idy),0,ierr)
              do iy=0,ny
              do ix=0,nx
                  wksp(ix,iy,idr)=wksp(ix,iy,idr)-alpha*wksp(ix,iy,idy)
              enddo
              enddo
          endif

          ctmp=dotProd(neqn,level0,ierr,wksp(0,0,idr),wksp(0,0,idr))
          rm=dsqrt(dreal(ctmp))/r0

          if(level.ge.1 .and. (m.eq.1.or.mod(m,10).eq.0)) then
                  print'(" rel_resid(",i4,")=",1pe8.2)',m,rm
          endif
          if(rm.le.tol)then
              if(level.ge.1) then
                  print'(" Rel_Resid(",i4,")=",1pe8.2)',m,rm
                  print*,"PCGNR converges in iteration=",m
              endif
              itCGNR=m
              return
          endif

          do iy=0,ny
          do ix=0,nx
              wksp(ix,iy,idz)=wksp(ix,iy,idr)
          enddo
          enddo
          call ilu9Subst(2,nx,ny,nbw,level0,ierr,ILUT,wksp(0,0,idz))
          call ilu9Subst(1,nx,ny,nbw,level0,ierr,ILU,wksp(0,0,idz))

          ctmp=dotProd(neqn,level0,ierr,wksp(0,0,idz),wksp(0,0,idr))
          rho1=dreal(ctmp)

          beta=rho1/rho0

          do iy=0,ny
          do ix=0,nx
              wksp(ix,iy,idp)= wksp(ix,iy,idz)+beta*wksp(ix,iy,idp)
          enddo
          enddo

          rho0=rho1

!--------------------------------
      enddo
!--------------------------------

      itCGNR=m
      return
      end


!=======================================================================
      subroutine ilu9(nx,ny,nbw,level,ierr,A,ILU)
!=======================================================================
      implicit none
      integer nx,ny,nbw,level,ierr
      complex*16 A(-nbw:nbw,0:nx,0:ny),ILU(-nbw:nbw,0:nx,0:ny)

      integer i,ix,iy
      complex*16 pivot,mul

      integer identity
      identity=0

!-------------------------
! Print User Information
!-------------------------

      if(level.ge.1) print'("ILU9: nx=",i4," ny=",i4)',nx,ny

      if(nbw.ne.4)then
          print*,"ERROR: ilu9.f: Wrong nbw=",nbw
          ierr=1
          return
      endif


!--------------------------------
! First, copy ILU <--- A
!--------------------------------

      do iy=0,ny
      do ix=0,nx
      do i =-nbw,nbw
          ILU(i,ix,iy)=A(i,ix,iy)
      enddo
      enddo
      enddo


!--------------------------------
! Now, ILU begins
!--------------------------------

      do iy=0,ny
      do ix=0,nx

          pivot=ILU(0,ix,iy)

          if(ix.ne.nx)then
                mul=ILU(-1,ix+1,iy)/pivot
                ILU(-1,ix+1,iy)=mul
                ILU( 0,ix+1,iy)=ILU( 0,ix+1,iy)-mul*ILU(1,ix,iy)
              if(iy.ne.ny)then
                ILU( 2,ix+1,iy)=ILU( 2,ix+1,iy)-mul*ILU(3,ix,iy)
                ILU( 3,ix+1,iy)=ILU( 3,ix+1,iy)-mul*ILU(4,ix,iy)
              endif
          endif

          if(iy.ne.ny)then

              if(ix.ne.0)then
                mul=ILU(-2,ix-1,iy+1)/pivot
                ILU(-2,ix-1,iy+1)=mul
                ILU( 0,ix-1,iy+1)=ILU( 0,ix-1,iy+1)-mul*ILU(2,ix,iy)
              if(ix.ne.nx)then
                ILU( 1,ix-1,iy+1)=ILU( 1,ix-1,iy+1)-mul*ILU(3,ix,iy)
              endif
              endif

                mul=ILU(-3,ix,iy+1)/pivot
                ILU(-3,ix,iy+1)=mul
              if(ix.ne.nx)then
                ILU(-2,ix,iy+1)=ILU(-2,ix,iy+1)-mul*ILU(1,ix,iy)
              endif
              if(ix.ne.0)then
                ILU(-1,ix,iy+1)=ILU(-1,ix,iy+1)-mul*ILU(2,ix,iy)
              endif
                ILU( 0,ix,iy+1)=ILU( 0,ix,iy+1)-mul*ILU(3,ix,iy)
              if(ix.ne.nx)then
                ILU( 1,ix,iy+1)=ILU( 1,ix,iy+1)-mul*ILU(4,ix,iy)
              endif

              if(ix.ne.nx)then
                mul=ILU(-4,ix+1,iy+1)/pivot
                ILU(-4,ix+1,iy+1)=mul
                ILU(-3,ix+1,iy+1)=ILU(-3,ix+1,iy+1)-mul*ILU(1,ix,iy)
              if(ix.ne.0)then
                ILU(-1,ix+1,iy+1)=ILU(-1,ix+1,iy+1)-mul*ILU(3,ix,iy)
              endif
                ILU( 0,ix+1,iy+1)=ILU( 0,ix+1,iy+1)-mul*ILU(4,ix,iy)
              endif
          endif

      enddo
      enddo

!===============================================
      if(identity.eq.1)then
      print*, "ilu9.f: Temp-Setting: Temp-Setting ILU=I"
      do iy=0,ny
      do ix=0,nx
          do i =-nbw,nbw
              ILU(i,ix,iy)=0.d0
          enddo
          ILU(0,ix,iy)=1.d0
      enddo
      enddo
      endif
!===============================================

      return
      end


!=======================================================================
      subroutine ilu9Subst(imode,nx,ny,nbw,level,ierr,ILU,x)
!=======================================================================
      implicit none
      integer imode,nx,ny,nbw,level,ierr
      complex*16 ILU(-nbw:nbw,0:nx,0:ny),x(0:nx,0:ny)

      integer ix,iy

! imode=1: Diag(L) = 1.d0
!       2: Diag(U) = 1.d0

!-------------------------
! Print User Information
!-------------------------

      if(level.ge.1) print'("ILU9SUBST: nx=",i4," ny=",i4)',nx,ny

      if(nbw.ne.4)then
          print*,"ERROR: ilu9Subst.f: Wrong nbw=",nbw
          ierr=1
          return
      endif


!---------------------------------
! forward ILU elimination: L^{-1}
!---------------------------------

      iy=0
          ix=0
              if(imode.eq.2) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          do ix=1,nx
              x(ix,iy)=x(ix,iy)-ILU(-1,ix,iy)*x(ix-1,iy)
              if(imode.eq.2) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          enddo

      do iy=1,ny
          ix=0
              x(ix,iy)=x(ix,iy)-ILU(-3,ix,iy)*x(ix,iy-1)
     &                         -ILU(-2,ix,iy)*x(ix+1,iy-1)
              if(imode.eq.2) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          do ix=1,nx-1
              x(ix,iy)=x(ix,iy)-ILU(-4,ix,iy)*x(ix-1,iy-1)
     &                         -ILU(-3,ix,iy)*x(ix,iy-1)
     &                         -ILU(-2,ix,iy)*x(ix+1,iy-1)
     &                         -ILU(-1,ix,iy)*x(ix-1,iy)
              if(imode.eq.2) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          enddo
          ix=nx
              x(ix,iy)=x(ix,iy)-ILU(-4,ix,iy)*x(ix-1,iy-1)
     &                         -ILU(-3,ix,iy)*x(ix,iy-1)
     &                         -ILU(-1,ix,iy)*x(ix-1,iy)
              if(imode.eq.2) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
      enddo


!-----------------------------------
! backward ILU substitution: U^{-1}
!-----------------------------------

      iy=ny
          ix=nx
              if(imode.eq.1) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          do ix=nx-1,0,-1
              x(ix,iy)=x(ix,iy)-ILU(1,ix,iy)*x(ix+1,iy)
              if(imode.eq.1) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          enddo

      do iy=ny-1,0,-1
          ix=nx
              x(ix,iy)=x(ix,iy)-ILU(2,ix,iy)*x(ix-1,iy+1)
     &                         -ILU(3,ix,iy)*x(ix,iy+1)
              if(imode.eq.1) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          do ix=nx-1,1,-1
              x(ix,iy)=x(ix,iy)-ILU(1,ix,iy)*x(ix+1,iy)
     &                         -ILU(2,ix,iy)*x(ix-1,iy+1)
     &                         -ILU(3,ix,iy)*x(ix,iy+1)
     &                         -ILU(4,ix,iy)*x(ix+1,iy+1)
              if(imode.eq.1) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
          enddo
          ix=0
              x(ix,iy)=x(ix,iy)-ILU(1,ix,iy)*x(ix+1,iy)
     &                         -ILU(3,ix,iy)*x(ix,iy+1)
     &                         -ILU(4,ix,iy)*x(ix+1,iy+1)
              if(imode.eq.1) x(ix,iy)=x(ix,iy)/ILU(0,ix,iy)
      enddo

      return
      end


!=======================================================================
      subroutine transConjg(nx,ny,nbw,level,ierr,A,B)
!=======================================================================
      implicit none
      integer nx,ny,nbw,level,ierr
      complex*16 A(-nbw:nbw,0:nx,0:ny),B(-nbw:nbw,0:nx,0:ny)

      integer i,ix,iy
      if(level.ge.2) print'("TransConjg: nx=",i4," ny=",i4)',nx,ny

!------------------
! initial setting
!------------------

      do iy=0,ny
      do ix=0,nx
      do i =-nbw,nbw
          B(i,ix,iy)=(0.d0,0.d0)
      enddo
      enddo
      enddo

!-----------------------
      if(nbw.eq.2)then
!-----------------------
          print*,"Error: transConjg.f: wrong nbw=",nbw
          ierr=1
          return

!-----------------------
      else if(nbw.eq.4)then
!-----------------------

      do iy=0,ny
      do ix=0,nx
          if(iy.ne.0) then
              if(ix.ne.0)then
                  B(-4,ix,iy)=dconjg(A( 4,ix-1,iy-1))
              endif
                  B(-3,ix,iy)=dconjg(A( 3,ix,iy-1))
              if(ix.ne.nx)then
                  B(-2,ix,iy)=dconjg(A( 2,ix+1,iy-1))
              endif
          endif
              if(ix.ne.0)then
                  B(-1,ix,iy)=dconjg(A( 1,ix-1,iy))
              endif
                  B( 0,ix,iy)=dconjg(A( 0,ix,iy))
              if(ix.ne.nx)then
                  B( 1,ix,iy)=dconjg(A(-1,ix+1,iy))
              endif
          if(iy.ne.nx) then
              if(ix.ne.0)then
                  B( 2,ix,iy)=dconjg(A(-2,ix-1,iy+1))
              endif
                  B( 3,ix,iy)=dconjg(A(-3,ix,iy+1))
              if(ix.ne.nx)then
                  B( 4,ix,iy)=dconjg(A(-4,ix+1,iy+1))
              endif
          endif
      enddo
      enddo

!-----------------------
      endif
!-----------------------

      return
      end


!=======================================================================
      subroutine mtx9vec(nx,ny,a,x,y,level,ierr)
!=======================================================================
! matrix-vector multiplication for row-compact matrix
! obtained from the 9-point finite difference method for
! 2D partial differential equations (PDE).
! INPUT:
!   a:  the matrix
!   x:  the vector
! OUTPUT:
!   y = a x
!----------------------------------------------------------------------
      implicit none
      integer nx,ny,level,ierr
      complex*16 a(-4:4,0:nx,0:ny),x(0:nx,0:ny),y(0:nx,0:ny)

      integer i,j,nxm1,nym1,im,ip,jm,jp

      if (level.ge.4) then
         print'("MTX9VEC: matrix vector multiplication")'
         print*,"dim(A)=",(nx+1)*(ny+1)
      endif

      nxm1=nx-1
      nym1=ny-1

!-------------------------
!---- interior points
!-------------------------

      do j=1,nym1
         jm=j-1
         jp=j+1
      do i=1,nxm1
         im=i-1
         ip=i+1
         y(i,j)=a(-4,i,j)*x(im,jm) +a(-3,i,j)*x(i ,jm)
     &         +a(-2,i,j)*x(ip,jm) +a(-1,i,j)*x(im,j )
     &         +a( 0,i,j)*x(i ,j ) +a( 1,i,j)*x(ip,j )
     &         +a( 2,i,j)*x(im,jp) +a( 3,i,j)*x(i ,jp)
     &         +a( 4,i,j)*x(ip,jp)
      end do
      end do

!-------------------
!---- four corners
!-------------------

      i=0
      j=0
         ip=i+1
         jp=j+1
      y(i,j)=a( 0,i,j)*x(i ,j ) +a( 1,i,j)*x(ip,j )
     &      +a( 3,i,j)*x(i ,jp) +a( 4,i,j)*x(ip,jp)

      i=nx
      j=0
         im=i-1
         jp=j+1
      y(i,j)=a(-1,i,j)*x(im,j ) +a( 0,i,j)*x(i ,j )
     &      +a( 2,i,j)*x(im,jp) +a( 3,i,j)*x(i ,jp)

      i=0
      j=ny
         ip=i+1
         jm=j-1
      y(i,j)=a(-3,i,j)*x(i ,jm) +a(-2,i,j)*x(ip,jm)
     &      +a( 0,i,j)*x(i ,j ) +a( 1,i,j)*x(ip,j )

      i=nx
      j=ny
         im=i-1
         jm=j-1
      y(i,j)=a(-4,i,j)*x(im,jm) +a(-3,i,j)*x(i ,jm)
     &      +a(-1,i,j)*x(im,j ) +a( 0,i,j)*x(i ,j )

!-------------------
!---- four sides
!-------------------

      i=0
         ip=i+1
      do j=1,nym1
         jm=j-1
         jp=j+1
         y(i,j)=a(-3,i,j)*x(i ,jm) +a(-2,i,j)*x(ip,jm)
     &         +a( 0,i,j)*x(i ,j ) +a( 1,i,j)*x(ip,j )
     &         +a( 3,i,j)*x(i ,jp) +a( 4,i,j)*x(ip,jp)
      end do

      i=nx
         im=i-1
      do j=1,nym1
         jm=j-1
         jp=j+1
         y(i,j)=a(-4,i,j)*x(im,jm) +a(-3,i,j)*x(i ,jm)
     &         +a(-1,i,j)*x(im,j ) +a( 0,i,j)*x(i ,j )
     &         +a( 2,i,j)*x(im,jp) +a( 3,i,j)*x(i ,jp)
      end do

      j=0
         jp=j+1
      do i=1,nxm1
         im=i-1
         ip=i+1
         y(i,j)=a(-1,i,j)*x(im,j ) +a( 0,i,j)*x(i ,j )
     &         +a( 1,i,j)*x(ip,j ) +a( 2,i,j)*x(im,jp)
     &         +a( 3,i,j)*x(i ,jp) +a( 4,i,j)*x(ip,jp)
      end do

      j=ny
         jm=j-1
      do i=1,nxm1
         im=i-1
         ip=i+1
         y(i,j)=a(-4,i,j)*x(im,jm) +a(-3,i,j)*x(i ,jm)
     &         +a(-2,i,j)*x(ip,jm) +a(-1,i,j)*x(im,j )
     &         +a( 0,i,j)*x(i ,j ) +a( 1,i,j)*x(ip,j )

      end do

      return
      end


!=======================================================================
      complex*16 function dotProd(neqn,level,ierr,x,y)
!=======================================================================
      implicit none
      integer neqn,level,ierr
      complex*16 x(neqn),y(neqn)

      integer i

!-------------------------
      if(level.ge.2) print'("DotProd: neqn=",i8)',neqn

      dotProd=0.d0

      do i=1,neqn
          dotProd=dotProd+dconjg(x(i))*y(i)
      enddo

      return
      end