! *****************************************************************
! *****************************************************************

subroutine algorithm1(n,x,f,g,gsupn,fmin,epsa,epsr,iter,fcnt,gcnt,hcnt,lsyscnt,inform)

  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: n
  integer, intent(out) :: fcnt,gcnt,hcnt,inform,iter,lsyscnt
  real(kind=8), intent(in) :: epsa,epsr,fmin
  real(kind=8), intent(out) :: f,gsupn

  ! ARRAY ARGUMENTS
  real(kind=8), intent(inout) :: x(n)
  real(kind=8), intent(out) :: g(n)

  ! LOCAL SCALARS
  integer :: flag
  real(kind=8), parameter :: L0 = 1.0d-06
  real(kind=8) :: fnew,gsupn0,L

  ! LOCAL ARRAYS
  real(kind=8) :: H(n,n),xnew(n)

  inform = 0

  iter = 0
  fcnt = 0
  gcnt = 0
  hcnt = 0
  lsyscnt = 0

  L = L0

  call evalf(n,x,f,flag)
  fcnt = fcnt + 1

  if (flag .ne. 0 ) then
     inform = - 1
     return
  end if

10 continue

  call evalg(n,x,g,flag)
  gcnt = gcnt + 1

  if (flag .ne. 0 ) then
     inform = - 2
     return
  end if

  gsupn = maxval( abs( g(1:n) ) )

  if ( iter .eq. 0 ) then
     gsupn0 = gsupn
  end if

  write(*,*) 'iter = ',iter,' f = ',f,' gsupn = ',gsupn

  open(20,file='kss-tabline-interrupted.txt')
  write(20,9100) n,inform,f,gsupn,iter,fcnt,gcnt,hcnt,lsyscnt,0.0d0
  close(20)

  if ( f .le. fmin ) then
     inform = 1
     return
  end if

  if ( gsupn .le. epsa ) then
     inform = 2
     return
  end if

  if ( gsupn .le. epsr * gsupn0 ) then
     inform = 3
     return
  end if

  iter = iter + 1

  call evalh(n,x,H,flag)
  hcnt = hcnt + 1

  if (flag .ne. 0 ) then
     inform = - 3
     return
  end if

  call compnewiterate(n,x,f,g,H,xnew,fnew,L,fcnt,lsyscnt,inform)
  if ( inform .ne. 0 ) return

  f = fnew

  x(1:n) = xnew(1:n)

  go to  10

9100 format(1X,I6,1X,I2,1X,1P,D14.6,1X,1P,D7.1,5(1X,I6),0P,F8.2)

end subroutine algorithm1

! *****************************************************************
! *****************************************************************

subroutine compnewiterate(n,x,f,g,H,xnew,fnew,L,fcnt,lsyscnt,inform)

  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: n
  integer, intent(inout) :: fcnt,inform,lsyscnt
  real(kind=8), intent(in) :: f
  real(kind=8), intent(inout) :: L
  real(kind=8), intent(out) :: fnew

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in) :: g(n),H(n,n),x(n)
  real(kind=8), intent(out) :: xnew(n)

  ! PARAMETERS
  real(kind=8), parameter :: gamma = 1.0d-04, eta = 0.25d0, delta = 1.0d-16

  ! LOCAL SCALARS
  integer :: dsyevinfo,flag
  real(kind=8) :: alpha,ared,atmp,gts,geucn,mu,pred

  ! LOCAL ARRAYS
  real(kind=8) :: dwork(10*n),lambda(n),s(n),shiftedl(n),y(n)

  ! =============================================================
  ! Step 0: Compute the spectral decomposition of H
  ! =============================================================

  call dsyev('V','U',n,H,n,lambda,dwork,10*n,dsyevinfo)

  if ( dsyevinfo .ne. 0 ) then
     write(*,*) 'dsyevinfo = ',dsyevinfo
     stop
  end if

  ! =============================================================
  ! Step 1: Compute muminus and muplus
  ! =============================================================

  geucn = norm2( g(1:n) )

  ! This is muminus (setting mu = muminus defines Algorithm S1)
  ! mu = 0.5d0 * ( sqrt( lambda(1) ** 2 + 4.0d0 * L * geucn ) - lambda(1) )

  ! This is muplus (setting mu = muminus defines Algorithm S2)
  mu = max( 0.0d0, - lambda(1) ) + sqrt( L * geucn )

  ! =============================================================
  ! Step 2: Solve the linear system
  ! =============================================================

  y(1:n) = - matmul( transpose( H(1:n,1:n) ), g(1:n) )

  shiftedl(1:n) = lambda(1:n) + mu

  if ( minval( shiftedl(1:n) ) .le. 0.0d0 ) then
     write(*,*) 'Matrix ( Hk + muk I ) is not positive definite.'
     stop
  end if

  where( shiftedl(1:n) .ne. 0.0d0 )
     s(1:n) = y(1:n) / shiftedl(1:n)
  elsewhere
     s(1:n) = 0.0d0
  end where

  s(1:n) = matmul( H(1:n,1:n), s(1:n) )

  lsyscnt = lsyscnt + 1

!!$    res = max( 0.0d0, maxval( abs( y(1:n) ), shiftedl(1:n) .eq. 0.0d0 ) ) / &
!!$          max( 1.0d0, maxval( abs( y(1:n) ) ) )
!!$    write(*,*) 'Residual = ',res

  ! =============================================================
  ! Step 3: Line search
  ! =============================================================

  gts = dot_product( g(1:n), s(1:n) )

  alpha = 1.0d0

300 continue

  xnew(1:n) = x(1:n) + alpha * s(1:n)

  call evalf(n,xnew,fnew,flag)
  fcnt = fcnt + 1

  if (flag .ne. 0 ) then
     inform = - 1
     return
  end if

  if ( .not. fnew .le. f + alpha * gamma * gts ) then
     atmp = ( - gts * alpha ** 2 ) / ( 2.0d0 * ( fnew - f - alpha * gts ) )

     if ( 0.1d0 * alpha .le. atmp .and. atmp .le. 0.9d0 * alpha ) then
        alpha = atmp
     else
        alpha = 0.5d0 * alpha
     end if

     go to 300
  end if

  ! =============================================================
  ! Step 4: Update estimation of L
  ! =============================================================

  if ( alpha .lt. 1.0d0 ) then
     L = 2.0d0 * L

  else
     ared = f - fnew
     pred = - 0.5d0 * gts

     if ( ared .gt. eta * pred ) then
        L = max( 0.5d0 * L, delta )
     end if
  end if

end subroutine compnewiterate
