diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index fd5da15f..c22186c0 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -11,7 +11,7 @@ set(EXAMPLES "primes" "BFS" "word_counts" "tokens" "filter" "linefit"
   "fast_fourier_transform" "bucketed_dijkstra" "push_relabel_max_flow" 
   "knuth_morris_pratt" "huffman_tree" "decision_tree_c45" "karatsuba" "suffix_tree"
   "2d_linear_program" "box_kdtree" "radix_tree" "ray_trace" "hash_map" "oct_tree"
-  "3d_range" "rectangle_intersection")
+  "3d_range" "rectangle_intersection" "multi_BFS") 
 
 function(add_example NAME)
   add_executable(${NAME} ${NAME}.cpp)
@@ -22,3 +22,5 @@ foreach(example ${EXAMPLES})
   add_example(${example})
 endforeach()
 
+add_executable(schedulers schedulers.cpp)
+target_link_libraries(schedulers PRIVATE parlay)
diff --git a/examples/helper/ligra_light.h b/examples/helper/ligra_light.h
index ea39b76b..b71fab64 100644
--- a/examples/helper/ligra_light.h
+++ b/examples/helper/ligra_light.h
@@ -47,6 +47,13 @@ struct vertex_subset {
     if (is_sparse) return sparse;
     else return parlay::pack_index<vertex>(dense);
   }
+  template <typename F>
+  void apply(const F& f) {
+    if (is_sparse)
+      parlay::parallel_for(0, sparse.size(), [&] (long i) {f(sparse[i]);});
+    else 
+      parlay::parallel_for(0, dense.size(), [&] (long i) {if (dense[i]) f(i);});
+  }
 };
 
 template <typename V>
diff --git a/examples/multi_BFS.cpp b/examples/multi_BFS.cpp
new file mode 100644
index 00000000..fc51fd96
--- /dev/null
+++ b/examples/multi_BFS.cpp
@@ -0,0 +1,46 @@
+#include <iostream>
+#include <string>
+
+#include <parlay/primitives.h>
+#include <parlay/sequence.h>
+#include <parlay/internal/get_time.h>
+
+#include "multi_BFS.h"
+#include "helper/graph_utils.h"
+
+// **************************************************************
+// Driver
+// **************************************************************
+using vertex = int;
+using nested_seq = parlay::sequence<parlay::sequence<vertex>>;
+using graph = nested_seq;
+using utils = graph_utils<vertex>;
+
+int main(int argc, char* argv[]) {
+  auto usage = "Usage: BFS <n> || BFS <filename>";
+  if (argc != 2) std::cout << usage << std::endl;
+  else {
+    long n = 0;
+    graph G;
+    try { n = std::stol(argv[1]); }
+    catch (...) {}
+    if (n == 0) {
+      G = utils::read_symmetric_graph_from_file(argv[1]);
+      n = G.size();
+    } else {
+      G = utils::rmat_graph(n, 20*n);
+    }
+    utils::print_graph_stats(G);
+    parlay::sequence<node_info> result;
+    parlay::internal::timer t("Time");
+    for (int i=0; i < 3; i++) {
+      result = multi_BFS(1, G);
+      t.next("BFS");
+    }
+
+    for (int i=0; i < 5; i++)
+      std::cout << std::hex << result[i].visited.load() << ", " << result[i].visited_prev.load() << ", " << (int) result[i].d.load() << std::endl;
+    //long visited = parlay::reduce(parlay::map(result, parlay::size_of()));
+    //std::cout << "num vertices visited: " << visited << std::endl;
+  }
+}
diff --git a/examples/multi_BFS.h b/examples/multi_BFS.h
new file mode 100644
index 00000000..21c5fa92
--- /dev/null
+++ b/examples/multi_BFS.h
@@ -0,0 +1,77 @@
+#include <atomic>
+#include <utility>
+
+#include <parlay/primitives.h>
+#include <parlay/sequence.h>
+#include <parlay/io.h>
+
+#include "helper/ligra_light.h"
+
+struct node_info {
+  std::atomic<ulong> visited;
+  std::atomic<ulong> visited_prev;
+  std::atomic<char> d;
+  node_info() : visited(0), visited_prev(0), d(0) {}
+};
+  
+template <typename vertex, typename graph>
+auto multi_BFS(vertex start, const graph& G) {
+  parlay::internal::timer t;
+  int round = 0;
+  auto sizes = parlay::map(G, parlay::size_of());
+  t.next("map");
+  
+  vertex i = parlay::max_element(sizes) - sizes.begin();
+  t.next("max");
+  auto ngh = G[i];
+  if (ngh.size() < 63) abort();
+  auto sub_sizes = parlay::map(ngh, [&] (vertex j) {return std::make_pair(G[j].size(), j);});
+  sub_sizes = parlay::sort(sub_sizes, [&] (auto a, auto b) {return a.first > b.first;});
+  
+  auto vertices = parlay::tabulate<vertex>(64, [&] (long j) {
+						 if (j == 0) return i;
+						 else return sub_sizes[j-1].second;});
+  t.next("ngh");
+  
+  parlay::sequence<node_info> nodes(G.size());
+  for (int i = 0; i < 64; i++) 
+    nodes[vertices[i]].visited = 1ul << i;
+  t.next("init");
+  
+  auto edge_f = [&] (vertex u, vertex v) -> bool {
+		  ulong u_visited = nodes[u].visited.load();
+		  ulong v_visited = nodes[v].visited_prev.load();
+		  if ((u_visited | v_visited) != v_visited) {
+		    nodes[v].visited_prev.fetch_or(u_visited);
+		    char old_d = nodes[v].d.load();
+		    return (old_d < round &&
+			    nodes[v].d.compare_exchange_strong(old_d, round));
+		  } else return false; };
+
+  auto cond_f = [&] (vertex v) {
+		  return !(nodes[v].visited.load() & 1);};
+
+  auto frontier_map = ligra::edge_map(G, G, edge_f, cond_f);
+  auto frontier = ligra::vertex_subset<vertex>();
+  frontier.add_vertices(vertices);
+  t.next("head");
+
+  long total = 0;
+  while (frontier.size() > 0) {
+    round++;
+    long m = frontier.size();
+    total += m;
+    std::cout << "frontier size: " << std::dec <<  m << ", " << total << std::endl;
+      frontier = frontier_map(frontier);
+      t.next("map");
+      frontier.apply([&] (vertex v) {
+	if ((nodes[v].visited_prev.load() & 1) == 0) {
+	  ulong tmp = nodes[v].visited_prev.load();
+	  nodes[v].visited_prev = nodes[v].visited.load();
+	  nodes[v].visited = tmp;
+	} else nodes[v].visited = nodes[v].visited_prev.load();});
+      t.next("update");
+  }
+
+  return nodes;
+}