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

subroutine algorithm31(n,x,f,g,gsupn,lammin,fmin,epsga,epsgr,epsha, &
     epshr,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) :: epsga,epsgr,epsha,epshr,fmin
  real(kind=8), intent(out) :: f,gsupn,lammin

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

  ! LOCAL SCALARS
  integer :: dsyevinfo,flag
  real(kind=8) :: fnew,gsupn0,hnorm0

  ! LOCAL ARRAYS
  real(kind=8) :: dwork(10*n),H(n,n),lambda(n),xnew(n)
  
  inform = 0

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

  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) ) )

  ! This should be commented in order to use the 2nd-order stopping criteria (part 1)
  if ( f .le. fmin ) then
     inform = 1
     return
  end if

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

  ! This should be commented in order to use the 2nd-order stopping criteria (part 2)
  if ( gsupn .le. epsga ) then
     inform = 2
     return
  end if

  ! This should be commented in order to use the 2nd-order stopping criteria (part 3)
  if ( gsupn .le. epsgr * gsupn0 ) then
     inform = 3
     return
  end if

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

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

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

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

  lammin = lambda(1)

  if ( iter .eq. 0 ) then
     hnorm0 = maxval( abs( lambda(1:n) ) )
  end if
  
  write(*,*) 'iter = ',iter,' f = ',f,' gsupn = ',gsupn,' lambda1 = ',lammin

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

  ! This should be UNcommented in order to use the 2nd-order stopping criteria

!!$  if ( f .le. fmin ) then
!!$     inform = 1
!!$     return
!!$  end if
!!$
!!$  if ( gsupn .le. epsga .and. lammin .ge. - epsha ) then
!!$     inform = 2
!!$     return
!!$  end if
!!$
!!$  if ( gsupn .le. epsgr * gsupn0 .and. lammin .ge. - epshr * hnorm0 ) then
!!$     inform = 3
!!$     return
!!$  end if

  iter = iter + 1

  call algorithm32(n,x,f,g,H,lambda,xnew,fnew,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,1X,1P,D8.1,5(1X,I6),0P,F8.2)

end subroutine algorithm31

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

subroutine algorithm32(n,x,f,g,H,lambda,xnew,fnew,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(out) :: fnew

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

  ! PARAMETERS
  logical, parameter :: verbose = .true.
  real(kind=8), parameter :: alpha = 1.0d-08, M = 1.0d+03

  ! LOCAL SCALARS
  integer :: flag
  real(kind=8) :: lambar,mu,mu1,mu2,quo,rho,sminnorm,snorm,t,res

  ! LOCAL ARRAYS
  real(kind=8) :: q1(n),s(n),smin(n),shiftedl(n),y(n)

  ! =============================================================
  ! Step 1: Check whether system (6) with mu=0 is compatible or not
  ! =============================================================

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

  lambar = max( 0.0d0, - lambda(1) )

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

  if ( verbose ) then
     write(*,*) 'At Step 1:'
     write(*,*) 'lambar = max( 0.0, - lambda_1 ) = ',lambar
     write(*,*) 'Shifted eigenvalues (first 5 components): ',shiftedl(1:min(5,n))
     write(*,*) 'RHS y = -Q^T g (first 5 components) = ',y(1:min(5,n))
  end if

  ! Compute residual
  res = max( 0.0d0, maxval( abs( y(1:n) ), shiftedl(1:n) .eq. 0.0d0 ) ) / &
       max( 1.0d0, maxval( abs( y(1:n) ) ) )
  if ( res .gt. 1.0d-08 ) then
     if ( verbose ) then
        write(*,*) 'The linear system is NOT compatible.'
        write(*,*) 'residual = ',res
     end if
     rho = 0.0d0
     go to 500
  end if

!!$    ! This criterion needs some relative tolerance in order to be useful in practice
!!$    if ( count( shiftedl(1:n) .eq. 0.0d0 .and. abs( y(1:n) ) .gt. 1.0d-08 ) .gt. 0 ) then
!!$       if ( verbose ) then
!!$          write(*,*) 'The linear system is NOT compatible.'
!!$       end if
!!$       rho = 0.0d0
!!$       go to 500
!!$    end if

  if ( verbose ) then
     write(*,*) 'The linear system is compatible.'
     write(*,*) 'residual = ',res
  end if

  ! =============================================================
  ! Step 2: Compute minimum norm solution to (6) with mu=0
  ! =============================================================

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

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

  sminnorm = norm2( smin(1:n) )

  lsyscnt = lsyscnt + 1

  if ( sminnorm .eq. 0.0d0 ) then
     if ( lambar .gt. 0.0d0 ) then
        rho = 2.0d0 * M
     else ! if ( lambar .eq. 0.0d0 ) then
        rho = 0.0d0
     end if
  else
     rho = lambar / ( 3.0d0 *  sminnorm )
  end if
  
  if ( verbose ) then
     write(*,*) 'At Step2:'
     write(*,*) 'The minimum norm solution is given by (first 5 components) = ', &
          smin(1:min(5,n)),' (norm = ',sminnorm,')'
     write(*,*) 'rho = ',rho,' (If smaller than or equal to M, we go to Step 4.'
  end if

  if ( rho .le. M ) then
     go to 400
  end if

  ! =============================================================
  ! Step 3
  ! =============================================================

  q1 = H(1:n,1)

  snorm = lambar / ( 3.0d0 * M )

  if ( verbose ) then
     write(*,*) 'At Step 3:'
     write(*,*) 'Eigenvector q_1 associated with lambda_1 (first 5 components) = ',q1(1:min(5,n))
  end if

  ! Step 3.1

310 continue

  t = sqrt( snorm ** 2 - sminnorm ** 2 )

  s(1:n) = smin(1:n) + t * q1(1:n)

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

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

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

  if ( verbose ) then
     write(*,*) 'Desired norm = ',snorm,' t = ',t,' s (first 5 components) = ', &
          s(1:min(5,n)),' actual norm = ',norm2( s(1:n) )
  end if

  if ( fnew .le. f - alpha * snorm ** 3 ) then
     if ( verbose ) then
        write(*,*) 'xplus = x + s gave sufficient decrease'
     end if
     return
  end if

  ! Step 3.2

  if ( snorm .ge. 2.0d0 * sminnorm ) then
     snorm = 0.5d0 * snorm
     go to 310
  end if

  ! =============================================================
  ! Step 4
  ! =============================================================

400 continue

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

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

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

  if ( verbose ) then
     write(*,*) 'At Step 4:'
  end if

  if ( fnew .le. f - alpha * sminnorm ** 3 ) then
     if ( verbose ) then
        write(*,*) 'xplus = x + smin gave sufficient decrease'
        return
     end if
     return
  end if

  if ( verbose ) then
     write(*,*) 'xplus = x + smin did NOT give sufficient decrease. ftrial = ',fnew
  end if

  ! =============================================================
  ! Step 5
  ! =============================================================

500 continue

  rho = max( 0.1d0, rho )

  ! Step 5.1

510 continue

  if ( verbose ) then
     write(*,*) 'At Step 5:'
     write(*,*) 'rho = ',rho
  end if

  mu      = 1.0d0
  s(1:n)  = matmul( H(1:n,1:n), y(1:n) / ( shiftedl(1:n) + mu ) )
  snorm   = norm2( s(1:n) )
  quo     = ( lambar + mu ) / ( 3.0d0 * snorm )
  lsyscnt = lsyscnt + 1

  if ( quo .lt. rho ) then
     mu1 = mu

511  continue
     if ( mu .ge. 1.0d+20 ) then
        write(*,*) 'Something went wrong in the bracketing phase: mu is too large.'
        stop
     end if

     mu      = 2.0d0 * mu
     s(1:n)  = matmul( H(1:n,1:n), y(1:n) / ( shiftedl(1:n) + mu ) )
     snorm   = norm2( s(1:n) )
     quo     = ( lambar + mu ) / ( 3.0d0 * snorm )
     lsyscnt = lsyscnt + 1

     if ( quo .lt. rho ) then
        go to 511
     else if ( quo .le. 100.0d0 * rho ) then
        if ( verbose ) then
           write(*,*) 'We found mu = ',mu,' and s (first 5 components) = ',s(1:min(5,n))
        end if
        go to 520
     else ! this is the case quo .gt. 100.0d0 * rho
        mu2 = mu
        go to 513
     end if

  else if ( quo .le. 100.0d0 * rho ) then
     go to 520

  else ! this is the case quo .gt. 100.0d0 * rho
     mu2 = mu

512  continue
     if ( mu .le. 1.0d-20 ) then
        write(*,*) 'Something went wrong in the bracketing phase: mu is too smal.'
        stop
     end if

     mu      = 0.5d0 * mu
     s(1:n)  = matmul( H(1:n,1:n), y(1:n) / ( shiftedl(1:n) + mu ) )
     snorm   = norm2( s(1:n) )
     quo     = ( lambar + mu ) / ( 3.0d0 * snorm )
     lsyscnt = lsyscnt + 1

     if ( quo .gt. 100.0d0 * rho ) then
        go to 512
     else if ( quo .ge. rho ) then
        if ( verbose ) then
           write(*,*) 'We found mu = ',mu,' and s (first 5 components) = ',s(1:min(5,n))
        end if
        go to 520
     else ! this is the case quo .lt. rho
        mu1 = mu
        go to 513
     end if
  end if

513 continue
  if ( mu1 .gt. mu2 .or. mu2 - mu1 .le. 1.0d-16 ) then
     write(*,*) 'Something went wrong in the bissection: too small interval.'
     stop
  end if

  mu      = 0.5d0 * ( mu1 + mu2 )
  s(1:n)  = matmul( H(1:n,1:n), y(1:n) / ( shiftedl(1:n) + mu ) )
  snorm   = norm2( s(1:n) )
  quo     = ( lambar + mu ) / ( 3.0d0 * snorm )
  lsyscnt = lsyscnt + 1

  if ( quo .lt. rho ) then
     mu1 = mu
     go to 513

  else if ( quo .le. 100.0d0 * rho ) then
     if ( verbose ) then
        write(*,*) 'We found mu = ',mu,' and s (first 5 components) = ',s(1:min(5,n))
     end if
     go to 520

  else ! this is the case quo .gt. 100.0d0 * rho
     mu2 = mu
     go to 513
  end if

  ! Step 5.2

520 continue

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

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

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

  if ( fnew .le. f - alpha * snorm ** 3 ) then
     if ( verbose ) then
        write(*,*) 'xplus = x + s gave sufficient decrease'
     end if
     return
  end if

  if ( verbose ) then
     write(*,*) 'xplus = x + s did NOT give sufficient decrease'
  end if

  ! Step 5.3

  if ( mu .lt. 0.1d0 ) then
     if ( verbose ) then
        write(*,*) 'Since mu < 0.1, rho will be increased.'
     end if

     rho = 10.0d0 * quo
     go to 510
  end if

  ! =============================================================
  ! Step 6
  ! =============================================================

  if ( verbose ) then
     write(*,*) 'At Step 6: '
  end if

600 continue

  mu = 2.0d0 * mu

  s(1:n) = matmul( H(1:n,1:n), y(1:n) / ( shiftedl(1:n) + mu ) )

  lsyscnt = lsyscnt + 1

  if ( verbose ) then
     write(*,*) 'mu = ',mu,' and s (first 5 components) = ',s(1:min(5,n))
  end if

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

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

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

  if ( fnew .le. f - alpha * norm2( s(1:n) ) ** 3 ) then
     if ( verbose ) then
        write(*,*) 'xplus = x + s gave sufficient decrease'
     end if
     return
  end if

  go to 600

end subroutine algorithm32

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