diff --git a/ADOL-C/test/Makefile.am b/ADOL-C/test/Makefile.am
index 43d40a88..1b02a1f6 100644
--- a/ADOL-C/test/Makefile.am
+++ b/ADOL-C/test/Makefile.am
@@ -10,17 +10,22 @@
 ## 
 ##############################################################################
 
-noinst_PROGRAMS         = powexam speelpenning fminmax
+noinst_PROGRAMS         = powexam speelpenning fminmax parallel
 
 speelpenning_SOURCES    = speelpenning.cpp
 
 powexam_SOURCES         = powexam.cpp
 
+parallel_SOURCES        = parallel.cpp
+
 fminmax_SOURCES		= fminmax.cpp
 
 powexam_LDADD           = ../lib${adolclib}.la
 speelpenning_LDADD = ../lib${adolclib}.la
 fminmax_LDADD		= ../lib${adolclib}.la
+parallel_LDADD		= ../lib${adolclib}.la
+
+parallel_CXXFLAGS				= -pthread
 
 AM_CFLAGS               = @ac_adolc_cflags@
 AM_CXXFLAGS             = @ac_adolc_cxxflags@
@@ -35,7 +40,7 @@ EXTRA_DIST = run_tests
 # on Cygwin
 DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir)/inc
 
-test:	powexam speelpenning 
+test:	powexam speelpenning parallel
 	chmod u+x $(srcdir)/run_tests
 	$(srcdir)/run_tests
 
diff --git a/ADOL-C/test/parallel.cpp b/ADOL-C/test/parallel.cpp
new file mode 100644
index 00000000..b384b0d1
--- /dev/null
+++ b/ADOL-C/test/parallel.cpp
@@ -0,0 +1,96 @@
+#include <adolc/adouble.h>            // use of active doubles
+#include <adolc/drivers/drivers.h>    // use of "Easy to Use" drivers
+// gradient(.) and hessian(.)
+#include <adolc/taping.h>             // use of taping
+
+#include <iostream>
+using namespace std;
+
+#include <cstdlib>
+#include <thread>
+
+void derive(const short my_tape, double init)
+{
+  const int n = 100;
+
+  double *xp = new double[n];
+  double  yp = 0.0;
+  adouble *x = new adouble[n];        
+  adouble  y = 1;
+
+  for(int i = 0; i < n; i++)
+  {
+    xp[i] = (i + init)/(2. + i);           // some initialization
+  }
+
+  trace_on(my_tape);
+  
+  for(int i = 0; i < n; i++) {
+    x[i] <<= xp[i];
+    y *= x[i];
+  }
+
+  y >>= yp;
+
+  delete[] x;
+  
+  trace_off();
+
+  double* g = new double[n];
+  
+  gradient(my_tape, n, xp, g);
+
+  double** H = new double*[n];
+  
+  for(int i = 0; i < n; i++)
+  {
+    H[i] = (double*) new double[i + 1];
+  }
+  
+  hessian(my_tape, n, xp, H);
+
+  double errh = 0.;
+
+  for(int i = 0; i < n; i++)
+  {
+    for(int j = 0; j < n; j++)
+    {
+      if (i>j)
+      {
+        errh += fabs(H[i][j] - g[i]/xp[j]);
+      }
+    } // end for
+  } // end for
+
+  std::cout << "Computed Hessian in tape " << my_tape
+            << ", error = "
+            << errh
+            << std::endl;
+
+  for(int i = 0; i < n; ++i)
+  {
+    delete[] H[i];
+  }
+  
+  delete[] H;
+  delete[] g;
+}
+
+
+int main()
+{
+  std::vector<std::thread> threads;
+
+  for(int i = 1; i <= 10; ++i)
+  {
+    threads.push_back(std::thread(derive, i, (double) i));
+  }
+
+  for(auto& thread: threads)
+  {
+    thread.join();
+  }
+    
+  return 0;
+}
+