      PROGRAM testJFNK
      IMPLICIT NONE

      INTEGER:: nsize
      DOUBLE PRECISION, ALLOCATABLE:: psi(:)
      DOUBLE PRECISION, ALLOCATABLE:: func(:)
      DOUBLE PRECISION, ALLOCATABLE:: psi0(:)
      DOUBLE PRECISION, ALLOCATABLE:: delta_psi(:)
      DOUBLE PRECISION, ALLOCATABLE:: b(:)
      DOUBLE PRECISION, ALLOCATABLE:: temp_psi(:)
      DOUBLE PRECISION, ALLOCATABLE:: temp_func(:)
      INTEGER:: jfnk_loopi
      INTEGER:: jfnk_loop_max
      DOUBLE PRECISION:: jfnk_epsilon
      DOUBLE PRECISION:: ddot

      INTEGER:: rsize
      INTEGER:: nwork
      DOUBLE PRECISION, ALLOCATABLE:: work(:)
      INTEGER:: irc(5)
      INTEGER:: icntl(8)
      DOUBLE PRECISION:: cntl(5)
      INTEGER:: info(3)
      DOUBLE PRECISION:: rinfo(2)
      DOUBLE PRECISION:: Jacob(5,5)

      INTEGER matvec, precondLeft, precondRight, dotProd
      PARAMETER (matvec=1, precondLeft=2, precondRight=3, dotProd=4)
      INTEGER revcom, colx, coly, colz, nbscal

      INTEGER:: i,j,k,l
      double precision:: zero
      double precision:: one
      

      Jacob(:,:)=0.D0
      do i=1,5,1
         Jacob(i,i)=-1.D0
      enddo
      do i=2,5,1
         Jacob(i,i-1)=2.D0
      enddo
      
      zero=0.D0
      one =1.D0
      jfnk_epsilon = 1.D-5
      nsize = 5
      jfnk_loop_max = 2

      allocate(psi(nsize))
      allocate(psi0(nsize))
      allocate(delta_psi(nsize))
      allocate(func(nsize))
      allocate(temp_psi(nsize))
      allocate(temp_func(nsize))
      psi0(:) = 0.01D0
      psi(:) = psi0(:)


      rsize = nsize
      nwork = rsize*rsize+rsize*(nsize+5)+6*nsize+1
      write(10,*) 'Unknown space size (nsize) is',nsize
      write(10,*) 'Work space size (nwork) is',nwork
      allocate(work(nwork))
      work(:) = 0.D0


      call init_DGMRES(icntl,cntl)
* Error message stream, default is 6
      icntl(1) = 8
* Warning message stream, default is 6
      icntl(2) = 8
* Convergence history steam, default is 0 (no display)
      icntl(3) = 8
* Controls the location of preconditioning, no default
* 0 means no preconditioning
      icntl(4) = 0
* Determines the orthogonalization scheme, default is 0 (MGS)
      icntl(5) = 1
* 0 means no initial guess is supplied, default is 0
      icntl(6) = 0
* Maximum number of iterations, no default
      icntl(7) = nsize*nsize
* Controls residual calculation, default is 0
      icntl(8) = 0

* Convergence tolerance, default is 1.D-5
      cntl(1) = 1.D-5
* Normalizing factor alpha, default is 0.D0
      cntl(2) = 1.D0
* Normalizing factor beta, default is 0.D0
      cntl(3) = 1.D0
* Normalizing factor alpha**p, default is 0.D0
      cntl(4) = 0.D0
* Normalizing factor beta**p, default is 0.D0
      cntl(5) = 0.D0



*****************************************************************************
* start jfnk loop
*****************************************************************************
      DO jfnk_loopi=1,jfnk_loop_max

      CALL setMyFunc(nsize,psi,func)
      work(nsize+1:2*nsize) = func(1:nsize)
      write(10,*) '*'
      write(10,*) '*'
      write(10,*) '*'
      write(10,*) 'Starting JFNK loop',jfnk_loopi
      write(10,*) 'psi:',psi
      write(10,*) 'func:',func
      write(8,*) '*'
      write(8,*) '*'
      write(8,*) '*'
      write(8,*) 'Starting JFNK loop',jfnk_loopi
      write(8,*) 'psi:',psi
      write(8,*) 'func:',func


10    CALL drive_DGMRES(nsize,nsize,rsize,nwork,work,irc,icntl,cntl,
     &  info,rinfo)

      revcom = irc(1)
      colx = irc(2)
      coly = irc(3)
      colz = irc(4)
      nbscal = irc(5)

      write(10,*) 'irc:',revcom,colx,coly,colz,nbscal


      IF (revcom.EQ.matvec)THEN
* perform the matrix vector product
*  work(colz) <- A * work(colx)
* which, for JFNK, becomes
*  work(colz) <- (F(psi+work(colx))-F(psi))/epsilon
      write(10,*) 'Performing matrix-vector product'
      write(10,*) 'psi:',psi
*      write(10,*) 'work:',work
*      write(10,*) 'work(colx):',work(colx)
      write(10,*) 'work(colx:colx+nsize-1):',work(colx:colx+nsize-1)
*      write(10,*) 'work(colz):',work(colz)
      temp_psi(:) = psi(:)+jfnk_epsilon*work(colx:colx+nsize-1)
      write(10,*) 'temp_psi:',temp_psi
      temp_func(:)=0.D0
      CALL setMyFunc(nsize,temp_psi,temp_func)
      temp_func(:) = (temp_func(:)-func(:))/jfnk_epsilon
      write(10,*) 'temp_func:',temp_func
      temp_func(:)=matmul(jacob,work(colx:colx+nsize-1))
      write(10,*) 'temp_func:',temp_func
      work(colz:colz+nsize-1) = temp_func(:)
      write(10,*) 'work(colz:colz+nsize-1):',work(colz:colz+nsize-1)
      GOTO 10


      ELSE IF (revcom.EQ.precondLeft) THEN
* perform the left preconditioning
*  work(colz) <- M_1^{-1} * work(colx)
*      CALL dcopy(nsize,work(colx),1,work(colz),1)
      write(6,*) 'ERROR: revcom asked for left preconditioning!'
      CALL exit(0)
*      GOTO 10


      ELSE IF (revcom.EQ.precondRight) THEN
* perform the right preconditioning
*  work(colz) <- M_2^{-1} * work(colx)
*      CALL dcopy(nsize,work(colx),1,work(colz),1)
      write(6,*) 'ERROR: revcom asked for right preconditioning!'
      CALL exit(0)
*     GOTO 10


      ELSE IF (revcom.EQ.dotProd) THEN
* perform the scalar product
*  work(colz) <- work(colx) work(coly)
* The statement to perform the dot products can be
* written in the compcat form.
      CALL dgemv('C',nsize,nbscal,ONE,work(colx),nsize,
     & work(coly),1,ZERO,work(colz),1)
* For the sake of simplicity we write it as a do-loop
*      DO i=0,nbscal-1
*         work(colz+i) = ddot(nsize,work(colx+i*nsize),1,
*     &    work(coly),1)
*      ENDDO
      GOTO 10

      ENDIF


      IF (info(1).EQ.0) THEN
         write(8,*) ' Normal exit'
         write(8,*) ' Convergence after ',info(2),' iterations'
         write(8,*) ' Backwards error - preconditioned system',
     &    rinfo(1)
         write(8,*) ' Backwards error - unpreconditioned system',
     &    rinfo(2)
         write(8,*) ' Solution : '
         DO j=1,nsize
            write(8,*) work(j)
         ENDDO
         write(8,*) ' Optimal size for workspace ',info(3)
      ELSE IF (info(1).EQ.-1) THEN
         write(8,*) ' Bad value of n'
      ELSE IF (info(1).EQ.-2) THEN
         write(8,*) ' Bad value of m'
      ELSE IF (info(1).EQ.-3) THEN
         write(8,*) ' Too small workspace'
         write(8,*) ' Minimal value should be ',info(2)
      ELSE IF (info(1).EQ.-4) THEN
         write(8,*) ' No convergence after ',icntl(7),' iterations'
      ELSE IF (info(1).EQ.-5) THEN
         write(8,*) ' Type of preconditioner not specified'
      ENDIF

      delta_psi(:) = work(1:nsize)
      psi(:) =  psi(:)-delta_psi(:)
      write(10,*) 'delta_psi:',delta_psi
      write(10,*) 'End of Newton loop:',jfnk_loopi,psi

      ENDDO
*****************************************************************************
* end jfnk loop
*****************************************************************************


      write(10,*) '*'
      write(10,*) '*'
      write(10,*) '*'
      write(10,*) 'Finished JFNK loops.'
      write(10,*) 'psi:',psi
      write(10,*) 'func:',func



      STOP
      END PROGRAM testJFNK



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


      SUBROUTINE setMyFunc(nsize,my_psi,my_func)
      IMPLICIT NONE
      INTEGER:: nsize
      DOUBLE PRECISION:: my_psi(nsize)
      DOUBLE PRECISION:: my_func(nsize)
      INTEGER:: i

      my_func(1) = 1.D0-my_psi(1)
      DO i=2,nsize
         my_func(i) = 2.D0*my_psi(i-1)-my_psi(i)
      ENDDO

      END SUBROUTINE setMyFunc
