fixed_point_linesearch Subroutine

    subroutine fixed_point_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)

    implicit none

    class(nlesolver_type),intent(inout) :: me
    real(wp),dimension(me%n),intent(in) :: xold      !! previous value of `x`
    real(wp),dimension(me%n),intent(in) :: p         !! search direction
    real(wp),dimension(me%n),intent(out) :: x        !! new `x`
    real(wp),intent(inout) :: f                      !! magnitude of `fvec`
    real(wp),dimension(me%m),intent(inout) :: fvec   !! function vector
    real(wp),dimension(:,:),intent(in),optional :: fjac !! jacobian matrix [dense]
    real(wp),dimension(:),intent(in),optional :: fjac_sparse !! jacobian matrix [sparse]

    integer :: i !! counter
    integer :: n_points !! number of points to compute
    real(wp),dimension(:),allocatable :: alphas_to_try !! set of `alpha` values to try
    real(wp),dimension(:),allocatable :: x_tmp !! temp `x`
    real(wp),dimension(:),allocatable :: fvec_tmp !! temp `fvec`
    real(wp) :: f_tmp !! temp `f`
    real(wp) :: step_size !! step size for `alpha`
    integer :: n !! number of steps to divide the interval

    ! 1 o-----------o
    ! 2 o-----o-----o
    ! 3 o---o---o---o

    n = me%n_intervals
    n_points = n + 1

    allocate(alphas_to_try(n_points))
    allocate(x_tmp(me%n))
    allocate(fvec_tmp(me%m))

    step_size = (me%alpha_max - me%alpha_min) / real(n,wp)

    ! compute the alphas:
    alphas_to_try(1) = me%alpha_min
    do i = 2, n
        alphas_to_try(i) = alphas_to_try(i-1) + step_size
    end do
    alphas_to_try(n_points) = me%alpha_max

    ! now compute the functions at these alphas:
    f = big
    do i = 1, n_points

        x_tmp = xold + p * alphas_to_try(i)

        ! evaluate the function at tthis point:
        call me%func(x_tmp,fvec_tmp)
        f_tmp = norm2(fvec_tmp)

        if (f_tmp<=f) then ! new best point
            x = x_tmp
            f = f_tmp
            fvec = fvec_tmp
        end if

    end do

    end subroutine fixed_point_linesearch