diff --git a/src/SortingAlgorithms.jl b/src/SortingAlgorithms.jl
index b9b0c8e..64ed633 100644
--- a/src/SortingAlgorithms.jl
+++ b/src/SortingAlgorithms.jl
@@ -9,12 +9,13 @@ using Base.Order
 import Base.Sort: sort!
 import DataStructures: heapify!, percolate_down!
 
-export HeapSort, TimSort, RadixSort, CombSort
+export HeapSort, TimSort, RadixSort, CombSort, BitonicSort
 
 struct HeapSortAlg  <: Algorithm end
 struct TimSortAlg   <: Algorithm end
 struct RadixSortAlg <: Algorithm end
 struct CombSortAlg  <: Algorithm end
+struct BitonicSortAlg  <: Algorithm end
 
 const HeapSort  = HeapSortAlg()
 const TimSort   = TimSortAlg()
@@ -46,6 +47,25 @@ Characteristics:
 """
 const CombSort  = CombSortAlg()
 
+"""
+    BitonicSort
+
+Indicates that a sorting function should use the bitonic mergesort algorithm.
+This algorithm performs a series of pre-determined comparisons, and tends to be very parallelizable.
+The algorithm effectively implements a sorting network based on merging bitonic sequences.
+
+Characteristics:
+ - *not stable* does not preserve the ordering of elements which
+   compare equal (e.g. "a" and "A" in a sort of letters which
+   ignores case).
+ - *in-place* in memory.
+ - *parallelizable* suitable for vectorization with SIMD instructions because
+it performs many independent comparisons.
+
+See wikipedia.org/wiki/Bitonic_sorter for more information.
+"""
+const BitonicSort  = BitonicSortAlg()
+
 
 ## Heap sort
 
@@ -609,9 +629,8 @@ function sort!(v::AbstractVector, lo::Int, hi::Int, ::CombSortAlg, o::Ordering)
     interval = (3 * (hi-lo+1)) >> 2
 
     while interval > 1
-        @inbounds for j in lo:hi-interval
-            a, b = v[j], v[j+interval]
-            v[j], v[j+interval] = lt(o, b, a) ? (b, a) : (a, b)
+        for j in lo:hi-interval
+            @inbounds comparator!(v, j, j+interval, o)
         end
         interval = (3 * interval) >> 2
     end
@@ -619,4 +638,53 @@ function sort!(v::AbstractVector, lo::Int, hi::Int, ::CombSortAlg, o::Ordering)
     return sort!(v, lo, hi, InsertionSort, o)
 end
 
+function sort!(v::AbstractVector, lo::Int, hi::Int, ::BitonicSortAlg, o::Ordering)
+    return bitonicsort!(view(v, lo:hi), o::Ordering)
+end
+
+function bitonicsort!(data, o::Ordering)
+    N = length(data)
+    for n in 1:intlog2(N)
+        bitonicfirststage!(data, Val(n), o::Ordering)
+        for m in n-1:-1:1
+            bitonicsecondstage!(data, Val(m), o::Ordering)
+        end
+    end
+    return data
+end
+
+function bitonicfirststage!(v, ::Val{Gap}, o::Ordering) where Gap
+    N = length(v)
+    gap = 1 << Gap
+    halfgap = 1 << (Gap - 1)
+    for i in 0:gap:N-1
+        firstj = max(0, i + gap - N)
+        for j in firstj:(halfgap-1)
+            ia = i + j
+            ib = i + gap - j - 1
+            @inbounds comparator!(v, ia + 1, ib + 1, o)
+        end
+    end
+end
+
+function bitonicsecondstage!(v, ::Val{Gap}, o::Ordering) where Gap
+    N = length(v)
+    gap = 1 << Gap
+    for i in 0:gap:N-1
+        lastj = min(N - 1, N - (i + gap >> 1 + 1))
+        for j in 0:lastj
+            ia = i + j
+            ib = i + j + gap >> 1
+            @inbounds comparator!(v, ia + 1, ib + 1, o)
+        end
+    end
+end
+
+intlog2(n) = (n > 1) ? 8sizeof(n-1)-leading_zeros(n-1) : 0
+
+Base.@propagate_inbounds function comparator!(v, i, j, o)
+    a, b = v[i], v[j]
+    v[i], v[j] = lt(o, b, a) ? (b, a) : (a, b)
+end
+
 end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index 7ea2780..0f4a558 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,7 +5,7 @@ using Random
 
 a = rand(1:10000, 1000)
 
-for alg in [TimSort, HeapSort, RadixSort, CombSort]
+for alg in [TimSort, HeapSort, RadixSort, CombSort, BitonicSort]
     b = sort(a, alg=alg)
     @test issorted(b)
     ix = sortperm(a, alg=alg)
@@ -85,7 +85,7 @@ for n in [0:10..., 100, 101, 1000, 1001]
         end
 
         # unstable algorithms
-        for alg in [HeapSort, CombSort]
+        for alg in [HeapSort, CombSort, BitonicSort]
             p = sortperm(v, alg=alg, order=ord)
             @test isperm(p)
             @test v[p] == si
@@ -99,7 +99,7 @@ for n in [0:10..., 100, 101, 1000, 1001]
 
     v = randn_with_nans(n,0.1)
     for ord in [Base.Order.Forward, Base.Order.Reverse],
-        alg in [TimSort, HeapSort, RadixSort, CombSort]
+        alg in [TimSort, HeapSort, RadixSort, CombSort, BitonicSort]
         # test float sorting with NaNs
         s = sort(v, alg=alg, order=ord)
         @test issorted(s, order=ord)