diff --git a/ex3_benchmarks/userset.h b/ex3_benchmarks/userset.h
new file mode 100644
index 0000000..a734e1f
--- /dev/null
+++ b/ex3_benchmarks/userset.h
@@ -0,0 +1,44 @@
+#ifndef USERSET_FILE
+#define USERSET_FILE
+#include <cmath>
+/**
+ * User function: f(@p x,@p y)
+ * @param[in] x		x-coordinate of discretization point
+ * @param[in] y		y-coordinate of discretization point
+ * @return  value for right hand side f(@p x,@p y)
+ */
+double FunctF(double const x, double const y);
+
+/**
+ * User function: u(@p x,@p y)
+ * @param[in] x		x-coordinate of discretization point
+ * @param[in] y		y-coordinate of discretization point
+ * @return  value for solution vector u(@p x,@p y)
+ */
+double FunctU(double const x, double const y);
+
+
+/**
+ * User function: f(@p x,@p y) = @f$ x^2 \sin(2.5\pi y)@f$.
+ * @param[in] x		x-coordinate of discretization point
+ * @param[in] y		y-coordinate of discretization point
+ * @return  value f(@p x,@p y)
+ */
+inline double fNice(double const x, double const y)
+{
+    return x * x * std::sin(2.5 * 3.14159 * y);
+}
+
+/**
+ * User function: f(@p x,@p y) = 0$.
+ * @param[in] x		x-coordinate of discretization point
+ * @param[in] y		y-coordinate of discretization point
+ * @return  value 0
+ */
+inline double f_zero(double const x, double const y)
+//double f_zero(double const  /*x*/, double const /*y*/)
+{
+    return 0.0 + 0.0*(x+y);
+}
+
+#endif
diff --git a/ex3_benchmarks/vdop.cpp b/ex3_benchmarks/vdop.cpp
new file mode 100644
index 0000000..55fac1a
--- /dev/null
+++ b/ex3_benchmarks/vdop.cpp
@@ -0,0 +1,84 @@
+#include "vdop.h"
+#include <cassert>               // assert()
+#include <cmath>
+#include <iostream>
+#include <vector>
+using namespace std;
+
+
+void vddiv(vector<double> &x, vector<double> const &y,
+           vector<double> const &z)
+{
+    assert( x.size() == y.size() && y.size() == z.size() );
+    size_t n = x.size();
+
+    for (size_t k = 0; k < n; ++k)
+    {
+        x[k] = y[k] / z[k];
+    }
+    return;
+}
+
+//******************************************************************************
+
+void vdaxpy(std::vector<double> &x, std::vector<double> const &y,
+            double alpha, std::vector<double> const &z )
+{
+    assert( x.size() == y.size() && y.size() == z.size() );
+    size_t n = x.size();
+
+    for (size_t k = 0; k < n; ++k)
+    {
+        x[k] = y[k] + alpha * z[k];
+    }
+    return;
+}
+//******************************************************************************
+
+double dscapr(std::vector<double> const &x, std::vector<double> const &y)
+{
+    assert( x.size() == y.size());
+    size_t n = x.size();
+
+    double    s = 0.0;
+    for (size_t k = 0; k < n; ++k)
+    {
+        s += x[k] * y[k];
+    }
+
+    return s;
+}
+
+//******************************************************************************
+void DebugVector(vector<double> const &v)
+{
+    cout << "\nVector  (nnode = " << v.size() << ")\n";
+    for (size_t j = 0; j < v.size(); ++j)
+    {
+        cout.setf(ios::right, ios::adjustfield);
+        cout << v[j] << "   ";
+    }
+    cout << endl;
+
+    return;
+}
+//******************************************************************************
+bool CompareVectors(std::vector<double> const &x, int const n, double const y[], double const eps)
+{
+    bool bn = (static_cast<int>(x.size()) == n);
+    if (!bn)
+    {
+        cout << "#########   Error: " << "number of elements" << endl;
+    }
+    //bool bv = equal(x.cbegin(),x.cend(),y);
+    bool bv = equal(x.cbegin(), x.cend(), y,
+                    [eps](double a, double b) -> bool
+    { return std::abs(a - b) < eps * (1.0 + 0.5 * (std::abs(a) + std::abs(a))); }
+                   );
+    if (!bv)
+    {
+        assert(static_cast<int>(x.size()) == n);
+        cout << "#########   Error: " << "values" << endl;
+    }
+    return bn && bv;
+}
diff --git a/ex3_benchmarks/vdop.h b/ex3_benchmarks/vdop.h
new file mode 100644
index 0000000..b2d0adb
--- /dev/null
+++ b/ex3_benchmarks/vdop.h
@@ -0,0 +1,58 @@
+#ifndef VDOP_FILE
+#define VDOP_FILE
+#include <vector>
+
+/** @brief  Element-wise vector divison x_k = y_k/z_k.
+ *
+ * @param[out] x  target vector
+ * @param[in]  y  source vector
+ * @param[in]  z  source vector
+ *
+ */
+void vddiv(std::vector<double> & x, std::vector<double> const& y,
+                                    std::vector<double> const& z);
+
+/** @brief  Element-wise daxpy operation x(k) = y(k) + alpha*z(k).
+ *
+ * @param[out] x  target vector
+ * @param[in]  y  source vector
+ * @param[in]  alpha  scalar
+ * @param[in]  z  source vector
+ *
+ */
+void vdaxpy(std::vector<double> & x, std::vector<double> const& y,
+                       double alpha, std::vector<double> const& z );
+
+
+/** @brief  Calculates the Euclidian inner product of two vectors.
+ *
+ * @param[in]  x vector
+ * @param[in]  y vector
+ * @return Euclidian inner product @f$\langle x,y \rangle@f$
+ *
+ */
+double dscapr(std::vector<double> const& x, std::vector<double> const& y);
+
+
+/**
+ * Print entries of a vector.
+ * @param[in] v	    vector values
+*/
+void DebugVector(std::vector<double> const &v);
+
+/** @brief  Compares an STL vector with POD vector.
+ *
+ * The accuracy criteria @f$ |x_k-y_k| < \varepsilon \left({1+0.5(|x_k|+|y_k|)}\right) @f$
+ * follows the book by
+ * <a href="https://www.springer.com/la/book/9783319446592">Stoyan/Baran</a>, p.8.
+ *
+ * @param[in]  x    STL vector
+ * @param[in]  n    length of POD vector
+ * @param[in]  y    POD vector
+ * @param[in]  eps  relative accuracy criteria (default := 0.0).
+ * @return true iff pairwise vector elements are relatively close to each other.
+ *
+ */
+bool CompareVectors(std::vector<double> const& x, int n, double const y[], double const eps=0.0);
+
+#endif