--- a/liboctave/MArray.cc
+++ b/liboctave/MArray.cc
@@ -2,6 +2,7 @@
 
 Copyright (C) 1993, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
               2007, 2008 John W. Eaton
+Copyright (C) 2009 VZLU Prague              
 
 This file is part of Octave.
 
@@ -53,6 +54,62 @@
   return 0;
 }
 
+template <class T>
+struct _idxadds_helper
+{
+  T *array;
+  T val;
+  _idxadds_helper (T *a, T v) : array (a), val (v) { }
+  void operator () (octave_idx_type i)
+    { array[i] += val; }
+};
+
+template <class T>
+struct _idxadda_helper
+{
+  T *array;
+  const T *vals;
+  _idxadda_helper (T *a, const T *v) : array (a), vals (v) { }
+  void operator () (octave_idx_type i)
+    { array[i] += *vals++; }
+};
+
+template <class T>
+void
+MArray<T>::idx_add (const idx_vector& idx, T val)
+{
+  octave_idx_type n = this->length ();
+  octave_idx_type ext = idx.extent (n);
+  if (ext > n)
+    {
+      this->resize (ext);
+      n = ext;
+    }
+
+  OCTAVE_QUIT;
+
+  octave_idx_type len = idx.length (n);
+  idx.loop (len, _idxadds_helper<T> (this->fortran_vec (), val));
+}
+
+template <class T>
+void
+MArray<T>::idx_add (const idx_vector& idx, const MArray<T>& vals)
+{
+  octave_idx_type n = this->length ();
+  octave_idx_type ext = idx.extent (n);
+  if (ext > n)
+    {
+      this->resize (ext);
+      n = ext;
+    }
+
+  OCTAVE_QUIT;
+
+  octave_idx_type len = std::min (idx.length (n), vals.length ());
+  idx.loop (len, _idxadda_helper<T> (this->fortran_vec (), vals.data ()));
+}
+
 // Element by element MArray by scalar ops.
 
 template <class T>