Merge Gael's experimental OpenMP parallelization support into Assign.h.

Eigen/Core

@@ -10,6 +10,13 @@
 #endif
 #endif
 
+#ifndef EIGEN_DONT_USE_OPENMP
+#ifdef _OPENMP
+#define EIGEN_USE_OPENMP
+#include <omp.h>
+#endif
+#endif
+
 #include <cstdlib>
 #include <cmath>
 #include <complex>

Eigen/src/Core/Assign.h

@@ -3,6 +3,7 @@
 //
 // Copyright (C) 2007 Michael Olbrich <michael.olbrich@gmx.net>
 // Copyright (C) 2006-2008 Benoit Jacob <jacob@math.jussieu.fr>
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -175,8 +176,28 @@ struct ei_operator_equals_impl<Derived, OtherDerived, false>
           unroll ? Derived::SizeAtCompileTime : Dynamic
           >::run(dst.derived(), src.derived());
       else
-        for(int i = 0; i < dst.size(); i++)
+      {
-          dst.coeffRef(i) = src.coeff(i);
+        #ifdef EIGEN_USE_OPENMPf
+        if(Derived::Flags & OtherDerived::Flags & LargeBit)
+        {
+          #ifdef __INTEL_COMPILER
+          #pragma omp parallel default(none) shared(other)
+          #else
+          #pragma omp parallel default(none)
+          #endif
+          {
+            #pragma omp for
+            for(int i = 0; i < dst.size(); i++)
+              dst.coeffRef(i) = src.coeff(i);
+          }
+        }
+        else
+        #endif // EIGEN_USE_OPENMP
+        {
+          for(int i = 0; i < dst.size(); i++)
+            dst.coeffRef(i) = src.coeff(i);
+        }
+      }
     }
     else // copying a matrix expression into a matrix
     {
@@ -192,18 +213,56 @@ struct ei_operator_equals_impl<Derived, OtherDerived, false>
       {
         if(Derived::ColsAtCompileTime == Dynamic || Derived::RowsAtCompileTime != Dynamic)
         {
-          // traverse in column-major order
+          #ifdef EIGEN_USE_OPENMP
-          for(int j = 0; j < dst.cols(); j++)
+          if(Derived::Flags & OtherDerived::Flags & LargeBit)
-            for(int i = 0; i < dst.rows(); i++)
+          {
-              dst.coeffRef(i, j) = src.coeff(i, j);
+            #ifdef __INTEL_COMPILER
+            #pragma omp parallel default(none) shared(other)
+            #else
+            #pragma omp parallel default(none)
+            #endif
+            {
+              #pragma omp for
+              for(int j = 0; j < dst.cols(); j++)
+                for(int i = 0; i < dst.rows(); i++)
+                  dst.coeffRef(i, j) = src.coeff(i, j);
+            }
+          }
+          else
+          #endif // EIGEN_USE_OPENMP
+          {
+            // traverse in column-major order
+            for(int j = 0; j < dst.cols(); j++)
+              for(int i = 0; i < dst.rows(); i++)
+                dst.coeffRef(i, j) = src.coeff(i, j);
+          }
         }
         else
         {
-          // traverse in row-major order
+          #ifdef EIGEN_USE_OPENMP
-          // in order to allow the compiler to unroll the inner loop
+          if(Derived::Flags & OtherDerived::Flags & LargeBit)
-          for(int i = 0; i < dst.rows(); i++)
+          {
-            for(int j = 0; j < dst.cols(); j++)
+            #ifdef __INTEL_COMPILER
-              dst.coeffRef(i, j) = src.coeff(i, j);
+            #pragma omp parallel default(none) shared(other)
+            #else
+            #pragma omp parallel default(none)
+            #endif
+            {
+              #pragma omp for
+              for(int i = 0; i < dst.rows(); i++)
+                for(int j = 0; j < dst.cols(); j++)
+                  dst.coeffRef(i, j) = src.coeff(i, j);
+            }
+          }
+          else
+          #endif // EIGEN_USE_OPENMP
+          {
+            // traverse in row-major order
+            // in order to allow the compiler to unroll the inner loop
+            for(int i = 0; i < dst.rows(); i++)
+              for(int j = 0; j < dst.cols(); j++)
+                dst.coeffRef(i, j) = src.coeff(i, j);
+          }
         }
       }
     }

Eigen/src/Core/Matrix.h

@@ -358,7 +358,8 @@ using Eigen::RowVector##SizeSuffix##TypeSuffix;
 EIGEN_USING_MATRIX_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 2) \
 EIGEN_USING_MATRIX_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 3) \
 EIGEN_USING_MATRIX_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 4) \
-EIGEN_USING_MATRIX_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, X)
+EIGEN_USING_MATRIX_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, X) \
+EIGEN_USING_MATRIX_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, XL)
 
 #define EIGEN_USING_MATRIX_TYPEDEFS \
 EIGEN_USING_MATRIX_TYPEDEFS_FOR_TYPE(i) \

bench/benchmarkX.cpp

@@ -5,18 +5,30 @@
 using namespace std;
 USING_PART_OF_NAMESPACE_EIGEN
 
+#ifndef MATTYPE
+#define MATTYPE MatrixXd
+#endif
+
+#ifndef MATSIZE
+#define MATSIZE 20
+#endif
+
+#ifndef REPEAT
+#define REPEAT 100000
+#endif
+
 int main(int argc, char *argv[])
 {
-	MatrixXd I = MatrixXd::identity(20,20);
+	MATTYPE I = MATTYPE::identity(MATSIZE,MATSIZE);
-	MatrixXd m(20,20);
+	MATTYPE m(MATSIZE,MATSIZE);
-	for(int i = 0; i < 20; i++) for(int j = 0; j < 20; j++)
+	for(int i = 0; i < MATSIZE; i++) for(int j = 0; j < MATSIZE; j++)
 	{
-		m(i,j) = 0.1 * (i+20*j);
+		m(i,j) = 0.1 * (i+MATSIZE*j)/MATSIZE;
 	}
-	for(int a = 0; a < 100000; a++)
+	for(int a = 0; a < REPEAT; a++)
 	{
 		m = I + 0.00005 * (m + m*m);
 	}
-	cout << m << endl;
+	cout << m(0,0) << endl;
 	return 0;
 }

bench/benchmarkXL.cpp

@@ -0,0 +1,34 @@
+// g++ -O3 -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX
+
+#include <Eigen/Core>
+
+using namespace std;
+USING_PART_OF_NAMESPACE_EIGEN
+
+#ifndef MATTYPE
+#define MATTYPE MatrixXLd
+#endif
+
+#ifndef MATSIZE
+#define MATSIZE 400
+#endif
+
+#ifndef REPEAT
+#define REPEAT 10000
+#endif
+
+int main(int argc, char *argv[])
+{
+	MATTYPE I = MATTYPE::identity(MATSIZE,MATSIZE);
+	MATTYPE m(MATSIZE,MATSIZE);
+	for(int i = 0; i < MATSIZE; i++) for(int j = 0; j < MATSIZE; j++)
+	{
+		m(i,j) = 0.1 * (i+MATSIZE*j)/MATSIZE;
+	}
+	for(int a = 0; a < REPEAT; a++)
+	{
+		m = I + 0.00005 * (m + m/4);
+	}
+	cout << m(0,0) << endl;
+	return 0;
+}