*************************************************************************
**       TEST PROGRAMME FOR THE various types of sparse solvers
*************************************************************************

      program validation
*
      implicit none
      integer lda, ldstrt, lwork
      parameter (lda = 1000, ldstrt = 60)
      parameter (lwork = ldstrt**2 + ldstrt*(lda+5) + 5*lda + 1)
*
      integer i, j, n, m
      integer revcom, colx, coly, colz, nbscal
      integer irc(5), icntl(8), info(3)
*
      integer matvec, precondLeft, precondRight, dotProd
      parameter (matvec=1, precondLeft=2, precondRight=3, dotProd=4)
*
      integer nout
*
      complex*16  a(lda,lda), work(lwork), xtrue(lda)
      real*8  cntl(5), rinfo(2)
*
      complex*16 ZERO, ONE
      parameter (ZERO = (0.0d0, 0.0d0), ONE = (1.0d0, 0.0d0))
*
*
***************************************************************
** Generate the test matrix a and set the right-hand side
** in positions (n+1) to 2n of the array work.
** The right-hand side is chosen such that the exact solution
** is the vector of all ones.
***************************************************************
*
      write(*,*) '***********************************************'
      write(*,*) 'This code is an example of use of'
      write(*,*) '2cd Order Richardson iteration'
      write(*,*) 'in double precision complex arithmetic'
      write(*,*) 'Results are written in output files'
      write(*,*) 'fort.40 and sol_zTestgmres.'
      write(*,*) '***********************************************'
      write(*,*)
      write(*,*) 'Matrix size ( < ', lda, ") ?"
      read(*,*) n
!     n = 19
      write(*,*) "Matrix size =", n
      if (n.gt.lda) then
        write(*,*) 'You are asking for a too large matrix'
        goto 100
      endif
*
! Solves A x = b, where A = grad^2 Poisson operator:
      do j = 1,n
        do i = 1,n
          a(i,j) = ZERO
        enddo
! desired solution for x :
        work(j) = ONE
      enddo
*
      do i = 1,n
        a(i,i) = (-2.d0, 0.d0)
      enddo
      do i = 1,n-1
        a(i,i+1) = (1.d0, 0.d0)
        a(i+1,i) = (1.d0, 0.d0)
      enddo
*
! compute b needed to give the desired solution for x:
      call ZGEMV('N',n,n,ONE,A,lda,work(1),1,ZERO,work(n+1),1)
*
      do j = 1,n
        work(j) = ZERO
      enddo
! or give an explicit b:
      do j = 1,n
         work(n+j) = ZERO
      enddo
      work(n+n/2+1) = -1

! I know the exact answer should be:
      do j = 1,n/2+1
         xtrue(j) = 0.5*j
      enddo
      do j = n/2+2, n
         xtrue(j) = xtrue(j-1)-0.5
      enddo
*
*********************************
** Choose the restart parameter
*********************************
*
!      write(*,*) 'Restart  <', ldstrt
!      read(*,*) m
*
*******************************************************
** Initialize the control parameters to default value
*******************************************************
*
*     call init_zgmres(icntl,cntl)
*
*************************
*c Tune some parameters
*************************
*
** Tolerance
*      cntl(1) = 1.d-10
** Save the convergence history in file fort.40
*      icntl(3) = 40
** No preconditioning
*      icntl(4) = 0
** ICGS orthogonalization
*      icntl(5) = 3
** Maximum number of iterations
*      icntl(7) = 1000
*! use recurrence formula to calculate residual:
*      icntl(8) = 0
*
*****************************************
** 1st or 2cd order Richardson iterative solution of A x = b:
*****************************************
*
      call riter2(n, A, lda, work(1), work(n+1), info, rinfo)
*
*******************************
* dump the solution on a file
*******************************
*
      nout = 11
      open(nout,FILE='sol_zTestgmres',STATUS='unknown')
!     write(nout,*) 'Restart : ', m
      write(nout,*) 'info(1) = ',info(1),'  info(2) = ',info(2)
      write(nout,*) 'rinfo(1) = ',rinfo(1),'  rinfo(2) = ',rinfo(2)
      write(nout,*) 'Optimal workspace = ', info(3)
      write(nout,*) 'Solution : '
      do j=1,n
        write(nout,*) work(j)
      enddo
      write(nout,*) '   '
*

      write(12,*) "# x Phi(x), Phi_true(x)"
      write(12,*) "0   0.0     0.0"
      do j=1,n
         write(12,*) j, real(work(j)), real(xtrue(j))
      enddo
      write(12,*) n+1," 0.0    0.0"

      write(*,*) "# iterations =", info(2), " err=", rinfo(2)
100    continue
*
      stop
      end

      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*16 A(lda,lda)
      complex*16 x(n), b(n)
      integer info(2)
      real*8 rinfo(2)

      complex*16 Ax(n), xp(n), xn(n), res(n)
      complex*16, parameter :: ONE = (1.0d0, 0.0d0)
      complex*16, 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

      subroutine riter2(n,A,lda,x,b, info, rinfo)
      ! solve A x = b using 2cd order Richardson iteration
      ! on input, x contains the initial guess.
      implicit none
      integer n, lda
      complex*16 A(lda,lda)
      complex*16 x(n), b(n)
      integer info(2)
      real*8 rinfo(2)

      complex*16 Ax(n), xp(n), xn(n), res(n)
      complex*16, parameter :: ONE = (1.0d0, 0.0d0)
      complex*16, parameter :: ZERO = (0.0d0, 0.0d0)
      real*8 dt, nu, err, pi, kmin, alpha, beta1, nudt_half
      integer iter, niter, maxiter

      dt=0.5*2.0
      nu=2*2*3.14159/n * 1.1
      maxiter = 1000
      pi = 2*abs(acos(0.0))

      ! Calculated optimal iteration parameters:
      kmin=2*pi/(n+1)
      alpha = 1-2*kmin
      beta1 = 1/4.+alpha/4+sqrt(alpha)/2
      
      nudt_half=(1-alpha)/(1+alpha)
      dt = sqrt(beta1*(1+nudt_half))
!     dt = (1+2*kmin)/(1+kmin)/2*(1-kmin+sqrt(1-2*kmin))
!     dt = sqrt(dt)
      nu=nudt_half*2/dt*0.33

! must be an error in the above calculation of the optimal parameters, because 
! I empirically find a better fit is obtained if I multiply this nu by 0.33...

      write(*,*) "alpha, beta1 =", alpha, beta1
      write(*,*) "dt, nu =", dt, nu
      write(*,*) "predicted error reduction / iteration =", 1/(1+nu*dt)
      write(*,*) "b=", b
      write(15,*) "# iteration #,   err"
      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 + (x-xp)*(1-nu*dt/2)/(1+nu*dt/2) + dt**2*res/(1+nu*dt/2)
         xp = x
         x = xn
         err=maxval(abs(res))
         write(15,*) iter, err
         if(err .le. 1.e-10) goto 20
      enddo

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

      write(*,*) "average convergence rate /iteration", err**(1./niter)
      return
      end subroutine riter2