cnvex Subroutine

    subroutine cnvex(n, a, h)

        !! compute the upper convex hull of the set (i,a(i)), i.e., the set of
        !! vertices (i_k,a(i_k)), k=1,2,...,m, such that the points (i,a(i)) lie
        !! below the straight lines passing through two consecutive vertices.
        !! the abscissae of the vertices of the convex hull equal the indices of
        !! the true  components of the logical output vector h.
        !! the used method requires o(nlog n) comparisons and is based on a
        !! divide-and-conquer technique. once the upper convex hull of two
        !! contiguous sets  (say, {(1,a(1)),(2,a(2)),...,(k,a(k))} and
        !! {(k,a(k)), (k+1,a(k+1)),...,(q,a(q))}) have been computed, then
        !! the upper convex hull of their union is provided by the subroutine
        !! cmerge. the program starts with sets made up by two consecutive
        !! points, which trivially constitute a convex hull, then obtains sets
        !! of 3,5,9... points,  up to  arrive at the entire set.
        !! the program uses the subroutine  cmerge; the subroutine cmerge uses
        !! the subroutines left, right and ctest. the latter tests the convexity
        !! of the angle formed by the points (i,a(i)), (j,a(j)), (k,a(k)) in the
        !! vertex (j,a(j)) up to within a given tolerance toler, where i<j<k.

        implicit none

        integer,intent(in) :: n
        real(wp) :: a(n)
        logical,intent(out) :: h(n)

        integer :: i, j, k, m, nj, jc

        do i = 1, n
            h(i) = .true.
        end do
        ! compute k such that n-2<=2**k<n-1
        k = int(log(n - 2.0_wp)/log(2.0_wp))
        if (2**(k + 1) <= (n - 2)) k = k + 1
        ! for each m=1,2,4,8,...,2**k, consider the nj pairs of consecutive
        ! sets made up by m+1 points having the common vertex
        ! (jc,a(jc)), where jc=m*(2*j+1)+1 and j=0,...,nj,
        ! nj=max(0,int((n-2-m)/(m+m))).
        ! compute the upper convex hull of their union by means of the
        ! subroutine cmerge
        m = 1
        do i = 0, k
            nj = max(0, int((n - 2 - m)/(m + m)))
            do j = 0, nj
                jc = (j + j + 1)*m + 1
                call cmerge(n, a, jc, m, h)
            end do
            m = m + m
        end do

    end subroutine cnvex