taskflow
diff --git a/examples/parallel_for.cpp b/examples/parallel_for.cpp
@@ -34,13 +34,44 @@ void for_each_index(int N) {
   executor.run(taskflow).get();
 }
 
+
+void for_each_index_nested() {
+	auto N = 10;
+
+	if (N < 0)
+	{
+		throw std::runtime_error("N must be non-negative");
+	}
+
+	int res;  // result
+
+	tf::Executor executor;
+	tf::Taskflow taskflow("nestedfor");
+
+	std::mutex mtx;
+
+	taskflow.for_each_index_nested(0,5,1,[&mtx](int i, tf::Subflow&sfi){
+		sfi.for_each_index(25,30,1,[i,&mtx](int j){
+			std::scoped_lock lock(mtx);
+			std::cout << "foreach_index_nested " << i << " "  << j << std::endl;
+		});
+	});
+
+	executor.run(taskflow).wait();
+	//taskflow.dump(std::cout);
+
+
+}
+
 // ----------------------------------------------------------------------------
 
 // Function: main
 int main() {
 
   for_each(100);
   for_each_index(100);
+  for_each_index_nested();
+
 
   return 0;
 }

diff --git a/taskflow/core/algorithm/for_each.hpp b/taskflow/core/algorithm/for_each.hpp
@@ -202,6 +202,100 @@ Task FlowBuilder::for_each_index(B&& beg, E&& end, S&& inc, C&& c){
   return task;
 }
 
+
+template <typename B, typename E, typename S, typename C>
+Task FlowBuilder::for_each_index_nested(B&& beg, E&& end, S&& inc, C&& c){
+
+  using I = stateful_index_t<B, E, S>;
+  using namespace std::string_literals;
+
+  Task task = emplace(
+  [b=std::forward<B>(beg), 
+   e=std::forward<E>(end), 
+   a=std::forward<S>(inc), 
+   c=std::forward<C>(c)] (Subflow& sf) mutable {
+
+    // fetch the iterator values
+    I beg = b;
+    I end = e;
+    I inc = a;
+
+    if(is_range_invalid(beg, end, inc)) {
+      TF_THROW("invalid range [", beg, ", ", end, ") with step size ", inc);
+    }
+
+    size_t chunk_size = 1;
+    size_t W = sf._executor.num_workers();
+    size_t N = distance(beg, end, inc);
+
+    // only myself - no need to spawn another graph
+    if(W <= 1 || N <= chunk_size) {
+      for(size_t x=0; x<N; x++, beg+=inc) {
+        c(beg, sf);
+      }
+      return;
+    }
+
+    if(N < W) {
+      W = N;
+    }
+
+    std::atomic<size_t> next(0);
+
+    for(size_t w=0; w<W; w++) {
+
+      sf.emplace([&next, beg, inc, N, chunk_size, W, &c] (Subflow&sfi) mutable {
+
+        size_t p1 = 2 * W * (chunk_size + 1);
+        double p2 = 0.5 / static_cast<double>(W);
+        size_t s0 = next.load(std::memory_order_relaxed);
+
+        while(s0 < N) {
+
+          size_t r = N - s0;
+
+          // find-grained
+          if(r < p1) {
+            while(1) { 
+              s0 = next.fetch_add(chunk_size, std::memory_order_relaxed);
+              if(s0 >= N) {
+                return;
+              }
+              size_t e0 = (chunk_size <= (N - s0)) ? s0 + chunk_size : N;
+              auto s = static_cast<I>(s0) * inc + beg;
+              for(size_t x=s0; x<e0; x++, s+=inc) {
+                c(s, sfi);
+              }
+            }
+            break;
+          }
+          // coarse-grained
+          else {
+            size_t q = static_cast<size_t>(p2 * r);
+            if(q < chunk_size) {
+              q = chunk_size;
+            }
+            size_t e0 = (q <= r) ? s0 + q : N;
+            if(next.compare_exchange_strong(s0, e0, std::memory_order_acquire,
+                                                    std::memory_order_relaxed)) {
+              auto s = static_cast<I>(s0) * inc + beg;
+              for(size_t x=s0; x<e0; x++, s+= inc) {
+                c(s,sfi);
+              }
+              s0 = next.load(std::memory_order_relaxed); 
+            }
+          }
+        } 
+      }).name("pfg_"s + std::to_string(w));
+
+    }
+
+    sf.join();
+  });  
+
+  return task;
+}
+
 }  // end of namespace tf -----------------------------------------------------
 
 

diff --git a/taskflow/core/flow_builder.hpp b/taskflow/core/flow_builder.hpp
@@ -369,6 +369,52 @@ class FlowBuilder {
     template <typename B, typename E, typename S, typename C>
     Task for_each_index(B&& first, E&& last, S&& step, C&& callable);
 
+    /**
+    @brief constructs an index-based parallel-for task and passes a new Subflow to each subtask
+
+    @tparam B beginning index type (must be integral)
+    @tparam E ending index type (must be integral)
+    @tparam S step type (must be integral)
+    @tparam C callable type
+
+    @param first index of the beginning (inclusive)
+    @param last index of the end (exclusive)
+    @param step step size 
+    @param callable a callable object to apply to each valid index
+
+    @return a tf::Task handle
+
+    The task spawns a subflow that applies the callable object to each index in the range <tt>[first, last)</tt> with the step size.
+
+    This method is equivalent to the parallel execution of the following loop:
+
+    @code{.cpp}
+    // case 1: step size is positive
+    for(auto i=first; i<last; i+=step) {
+      taskflow.emplace([i](tf::Subflow&sf)
+      {
+        callable(sf,i);
+      })
+    }
+
+    // case 2: step size is negative
+    for(auto i=first, i>last; i+=step) {
+      taskflow.emplace([i](tf::Subflow&sf)
+      {
+        callable(sf,i);
+      })
+    }
+    @endcode
+
+    Arguments are templated to enable stateful passing using std::reference_wrapper.
+    The callable needs to take two arguments the first of which is a non const
+    reference to a tf::Subflow and the second is a single argument of the integral index type.
+
+    Please refer to @ref ParallelIterations for details.
+    */
+    template <typename B, typename E, typename S, typename C>
+    Task for_each_index_nested(B&& beg, E&& end, S&& inc, C&& c);
+
     // ------------------------------------------------------------------------
     // reduction
     // ------------------------------------------------------------------------

diff --git a/unittests/algorithm.cpp b/unittests/algorithm.cpp
@@ -288,42 +288,69 @@ void stateful_for_each(unsigned W, TYPE) {
   std::vector<int> vec;
   std::atomic<int> counter {0};
 
-  for(size_t n = 0; n <= 150; n++) {
+  for(size_t third = 1; third <= 50; third++) {
+
+    size_t n = third * 3;
+
     for(size_t c=0; c<=17; c++) {
 
-    std::vector<int>::iterator beg, end;
-    size_t ibeg = 0, iend = 0;
-    size_t half = n/2;
+    using ILimit = std::pair<size_t,size_t>;
+    using Limit = std::pair<std::vector<int>::iterator,std::vector<int>::iterator>;
+    std::array<Limit,3> limits;
+    std::array<ILimit,3> ilimits;
+
 
     taskflow.clear();
+
 
     auto init = taskflow.emplace([&](){ 
       vec.resize(n);
       std::fill_n(vec.begin(), vec.size(), -1);
 
-      beg = vec.begin();
-      end = beg + half;
-
-      ibeg = half;
-      iend = n;
-
+      auto initL = [&](Limit & l, ILimit & il, size_t offset){
+        l.first = vec.begin()+offset;
+        l.second = vec.begin()+offset+third;
+
+        il.first = l.first-vec.begin();
+        il.second = l.second-vec.begin();
+      };
+      initL(limits[0],ilimits[0],0);
+      initL(limits[1],ilimits[1],third);
+      initL(limits[2],ilimits[2],2*third);
       counter = 0;
     });
+
 
-    tf::Task pf1, pf2;
+    tf::Task pf1, pf2, pf3;
 
     pf1 = taskflow.for_each(
-      std::ref(beg), std::ref(end), [&](int& i){
+      std::ref(limits[0].first), std::ref(limits[0].second), [&](int& i){
       counter++;
       i = 8;
     });
 
     pf2 = taskflow.for_each_index(
-      std::ref(ibeg), std::ref(iend), size_t{1}, [&] (size_t i) {
+      std::ref(ilimits[1].first), std::ref(ilimits[1].second), size_t{1}, [&] (size_t i) {
         counter++;
         vec[i] = -8;
     });
 
+    #if true
+      pf3 = taskflow.for_each_index_nested(
+        std::ref(ilimits[2].first), std::ref(ilimits[2].second), size_t{1}, [&] (size_t i, tf::Subflow&sf) {
+          sf.silent_async([&](){
+            counter++;
+            vec[i] = -16;
+          });
+      });
+    #else
+      pf3 = taskflow.for_each_index(
+        std::ref(ilimits[2].first), std::ref(ilimits[2].second), size_t{1}, [&] (size_t i) {
+          counter++;
+          vec[i] = -16;
+      });
+    #endif
+
     //switch (type) {
 
     //  case GUIDED:
@@ -369,20 +396,22 @@ void stateful_for_each(unsigned W, TYPE) {
     //  break;
     //}
 
-    init.precede(pf1, pf2);
+    init.precede(pf1,pf2,pf3);
 
     executor.run(taskflow).wait();
     REQUIRE(counter == n);
 
-    for(size_t i=0; i<half; ++i) {
-      REQUIRE(vec[i] == 8);
-      vec[i] = 0;
-    }
+    auto check = [&](size_t set, size_t val){
+      for(size_t i=ilimits[set].first; i<ilimits[set].second; ++i) {
+        REQUIRE(vec[i] == val);
+      }
+    };
+
+    check(0,8);
+    check(1,-8);
+    check(2,-16);
+
 
-    for(size_t i=half; i<n; ++i) {
-      REQUIRE(vec[i] == -8);
-      vec[i] = 0;
-    }
     }
   }
 }