diff --git a/src/include/arith.h b/src/include/arith.h
index c3d94b7..713823f 100644
--- a/src/include/arith.h
+++ b/src/include/arith.h
@@ -86,6 +86,11 @@ namespace MultiArrayTools
 	{
 	    return a1 + a2;
 	}
+
+    	static inline T& selfApply(T& a1, const T& a2)
+	{
+	    return a1 += a2;
+	}
     };
 
     template <typename T>
diff --git a/src/include/multi_array_operation.cc.h b/src/include/multi_array_operation.cc.h
index 4b1abed..86b77d4 100644
--- a/src/include/multi_array_operation.cc.h
+++ b/src/include/multi_array_operation.cc.h
@@ -97,6 +97,18 @@ namespace MultiArrayTools
 	return out;
     }
 
+    template <typename T, class OperationClass>
+    template <class AOp, class... Indices>
+    auto OperationBase<T,OperationClass>::oto(const AOp& aop,
+                                              const std::shared_ptr<Indices>&... inds) const
+	-> MultiArray<T,typename Indices::RangeType...>
+    {
+	MultiArray<T,typename Indices::RangeType...> out(inds->range()...,
+                                                         static_cast<T>(0));
+        AOp::selfApply(out(inds...), THIS());
+	return out;
+    }
+
     template <typename T, class OperationClass>
     template <typename R, class... Args>
     auto OperationBase<T,OperationClass>::a(const std::shared_ptr<function<R,T,typename Args::value_type...>>& ll,
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index a5a3dfe..143129a 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -69,6 +69,10 @@ namespace MultiArrayTools
 	auto to(const std::shared_ptr<Indices>&... inds) const
 	    -> MultiArray<T,typename Indices::RangeType...>;
 
+        template <class AOp, class... Indices>
+	auto oto(const AOp& aop, const std::shared_ptr<Indices>&... inds) const
+	    -> MultiArray<T,typename Indices::RangeType...>;
+
         template <typename R, class... Args> // Args = Operation Classes
         auto a(const std::shared_ptr<function<R,T,typename Args::value_type...>>& ll, const Args&... args) const
             -> Operation<R,function<R,T,typename Args::value_type...>,OperationClass, Args...>;