#include "complex_real.h"

!-----------------------------------------------------------------------
      subroutine riter1(n,A,lda,x,b, info, rinfo)
      ! solve A x = b using 1st order Richardson iteration
      ! on input, x contains the initial guess.
      implicit none
      integer n, lda
      COMPLEX_REAL A(lda,lda)
      COMPLEX_REAL x(n), b(n)
      integer info(2)
      real*8 rinfo(2)

      COMPLEX_REAL Ax(n), xp(n), xn(n), res(n)
      COMPLEX_REAL, parameter :: ONE = (1.0d0, 0.0d0)
      COMPLEX_REAL, parameter :: ZERO = (0.0d0, 0.0d0)
      real*8 dt, err
      integer iter, niter, maxiter

      dt=0.4
      maxiter = 100

      write(*,*) "b=", b

      do iter = 1, maxiter

!     perform the matrix vector product Ax = A * x
         call zgemv('N',n,n,ONE,A,lda,x,1,ZERO,Ax,1)
         res = Ax - b
!        write(*,*) "x =", x
!        write(*,*) "res =", res
         xn = x + dt*res
         x = xn
         err=maxval(abs(res))
         write(*,*) "err =", err
         if(err .le. 1.e-10) goto 20
      enddo

 20   niter = iter
      info(2)=niter
      rinfo(2)=err

      return
      end subroutine riter1