started optimizing the parallel code for small matrices and vectors

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs

@@ -67,15 +67,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha.IsZero())
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index]);
             }
             else if (alpha.IsOne())
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + x[index]);
             }
             else
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index]));
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + (alpha * x[index]));
             }
         }
 
@@ -95,15 +95,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha.IsZero())
             {
-                CommonParallel.For(0, x.Length, index => result[index] = Complex.Zero);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = Complex.Zero);
             }
             else if (alpha.IsOne())
             {
-                CommonParallel.For(0, x.Length, index => result[index] = x[index]);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = x[index]);
             }
             else
             {
-                CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; });
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => { result[index] = alpha * x[index]; });
             }
         }
 
@@ -172,7 +172,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] + y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] + y[i]);
         }
 
         /// <summary>
@@ -207,7 +207,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] - y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] - y[i]);
         }
 
         /// <summary>
@@ -242,7 +242,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] * y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] * y[i]);
         }
 
         /// <summary>
@@ -277,7 +277,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
+            CommonParallel.ElementFor(0, y.Length, y.Length, index => { result[index] = x[index] / y[index]; });
         }
 
         /// <summary>

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs

@@ -61,15 +61,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha.IsZero())
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index]);
             }
             else if (alpha.IsOne())
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + x[index]);
             }
             else
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index]));
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + (alpha * x[index]));
             }
         }
 
@@ -89,15 +89,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha.IsZero())
             {
-                CommonParallel.For(0, x.Length, index => result[index] = Complex32.Zero);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = Complex32.Zero);
             }
             else if (alpha.IsOne())
             {
-                CommonParallel.For(0, x.Length, index => result[index] = x[index]);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = x[index]);
             }
             else
             {
-                CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; });
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => { result[index] = alpha * x[index]; });
             }
         }
 
@@ -167,7 +167,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] + y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] + y[i]);
         }
 
         /// <summary>
@@ -202,7 +202,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] - y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] - y[i]);
         }
 
         /// <summary>
@@ -237,7 +237,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] * y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] * y[i]);
         }
 
         /// <summary>
@@ -272,7 +272,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
+            CommonParallel.ElementFor(0, y.Length, y.Length, index => { result[index] = x[index] / y[index]; });
         }
 
         /// <summary>

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs

@@ -61,15 +61,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha == 0.0)
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index]);
             }
             else if (alpha == 1.0)
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + x[index]);
             }
             else
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index]));
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + (alpha * x[index]));
             }
         }
 
@@ -89,15 +89,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha == 0.0)
             {
-                CommonParallel.For(0, x.Length, index => result[index] = 0.0);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = 0.0);
             }
             else if (alpha == 1.0)
             {
-                CommonParallel.For(0, x.Length, index => result[index] = x[index]);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = x[index]);
             }
             else
             {
-                CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; });
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => { result[index] = alpha * x[index]; });
             }
         }
 
@@ -167,7 +167,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; });
+            CommonParallel.ElementFor(0, y.Length, y.Length, index => { result[index] = x[index] + y[index]; });
         }
 
         /// <summary>
@@ -202,7 +202,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; });
+            CommonParallel.ElementFor(0, y.Length, y.Length, index => { result[index] = x[index] - y[index]; });
         }
 
         /// <summary>
@@ -237,7 +237,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; });
+            CommonParallel.ElementFor(0, y.Length, y.Length, index => { result[index] = x[index] * y[index]; });
         }
 
         /// <summary>
@@ -272,7 +272,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
+            CommonParallel.ElementFor(0, y.Length, y.Length, index => { result[index] = x[index] / y[index]; });
         }
 
         /// <summary>

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs

@@ -61,15 +61,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha == 0.0)
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index]);
             }
             else if (alpha == 1.0)
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + x[index]);
             }
             else
             {
-                CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index]));
+                CommonParallel.ElementFor(0, y.Length, y.Length, index => result[index] = y[index] + (alpha * x[index]));
             }
         }
 
@@ -89,15 +89,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha == 0.0)
             {
-                CommonParallel.For(0, x.Length, index => result[index] = 0.0f);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = 0.0f);
             }
             else if (alpha == 1.0)
             {
-                CommonParallel.For(0, x.Length, index => result[index] = x[index]);
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => result[index] = x[index]);
             }
             else
             {
-                CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; });
+                CommonParallel.ElementFor(0, x.Length, x.Length, index => { result[index] = alpha * x[index]; });
             }
         }
 
@@ -125,8 +125,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            float sum = 0;
-            CommonParallel.For(0, y.Length, index => sum += y[index] * x[index]);
+            var sum = 0.0f;
+
+            for (var index = 0; index < y.Length; index++)
+            {
+                sum += y[index] * x[index];
+            }
+
             return sum;
         }
 
@@ -162,7 +167,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] + y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] + y[i]);
         }
 
         /// <summary>
@@ -197,7 +202,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] - y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] - y[i]);
         }
 
         /// <summary>
@@ -232,7 +237,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, i => result[i] = x[i] * y[i]);
+            CommonParallel.ElementFor(0, y.Length, y.Length, i => result[i] = x[i] * y[i]);
         }
 
         /// <summary>
@@ -267,7 +272,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
+            CommonParallel.ElementFor(0, y.Length, y.Length, index => { result[index] = x[index] / y[index]; });
         }
         
         /// <summary>

src/Numerics/Control.cs

@@ -56,6 +56,11 @@ namespace MathNet.Numerics
         /// </summary>
         private static int _parallelizeOrder = 64;
 
+        /// <summary>
+        /// The default cutoff point for order size for the matrix multiply in linear algebra provider.
+        /// </summary>
+        private static int _parallelizeElements = 256;
+
         /// <summary>
         /// Initializes static members of the Control class.
         /// </summary>
@@ -161,5 +166,24 @@ namespace MathNet.Numerics
             }
         }
 
+        /// <summary>
+        /// Gets or sets the number of elements a vector or matrix must contain before we multiply threads.
+        /// </summary>
+        /// <value>Number of elements. Default is 256.</value>
+        public static int ParallelizeElements
+        {
+            get
+            {
+                return _parallelizeElements;
+            }
+
+            set
+            {
+                if (_parallelizeElements > 2)
+                {
+                    _parallelizeElements = value;
+                }
+            }
+        }
     }
 }

src/Numerics/Threading/CommonParallel.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2011 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -31,7 +31,6 @@
 namespace MathNet.Numerics.Threading
 {
     using System;
-    using System.Numerics;
 
 #if !SILVERLIGHT
     using System.Collections.Concurrent;
@@ -43,6 +42,26 @@ namespace MathNet.Numerics.Threading
     /// </summary>
     public static class CommonParallel
     {
+        /// <summary>
+        /// Executes a for loop in which iterations may run in parallel. 
+        /// </summary>
+        /// <param name="fromInclusive">The start index, inclusive.</param>
+        /// <param name="toExclusive">The end index, exclusive.</param>
+        /// <param name="numberOfElements">The number of elements that will be iterated over.</param>
+        /// <param name="body">The body to be invoked for each iteration.</param>
+        /// <exception cref="ArgumentNullException">The <paramref name="body"/> argument is <c>null</c>.</exception>
+        /// <exception cref="AggregateException">At least one invocation of the body threw an exception.</exception>
+        public static void ElementFor(int fromInclusive, int toExclusive, int numberOfElements, Action<int> body)
+        {
+            var parallel = true;
+            if (Control.DisableParallelization || Control.NumberOfParallelWorkerThreads < 2 || numberOfElements < Control.ParallelizeElements)
+            {
+                parallel = false;
+            }
+            
+            For(fromInclusive, toExclusive, body, parallel);
+        }
+
         /// <summary>
         /// Executes a for loop in which iterations may run in parallel. 
         /// </summary>
@@ -53,21 +72,37 @@ namespace MathNet.Numerics.Threading
         /// <exception cref="AggregateException">At least one invocation of the body threw an exception.</exception>
         public static void For(int fromInclusive, int toExclusive, Action<int> body)
         {
-#if SILVERLIGHT
-            Parallel.For(fromInclusive, toExclusive, body);
-#else
+            var parallel = true;
             if (Control.DisableParallelization || Control.NumberOfParallelWorkerThreads < 2)
             {
-                for (var index = fromInclusive; index < toExclusive; index++)
-                {
-                    body(index);
-                }
+                parallel = false;
             }
-            else
+
+            For(fromInclusive, toExclusive, body, parallel);
+        }
+
+        /// <summary>
+        /// Executes a for loop in which iterations may run in parallel. 
+        /// </summary>
+        /// <param name="fromInclusive">The start index, inclusive.</param>
+        /// <param name="toExclusive">The end index, exclusive.</param>
+        /// <param name="body">The body to be invoked for each iteration.</param>
+        /// <param name="parallel">Use multiple threads.</param>
+        /// <exception cref="ArgumentNullException">The <paramref name="body"/> argument is <c>null</c>.</exception>
+        /// <exception cref="AggregateException">At least one invocation of the body threw an exception.</exception>
+        public static void For(int fromInclusive, int toExclusive, Action<int> body, bool parallel)
+        {
+            if (parallel)
             {
+#if SILVERLIGHT
+                Parallel.For(fromInclusive, toExclusive, body);
+#else
                 Parallel.ForEach(
                     Partitioner.Create(fromInclusive, toExclusive),
-                    new ParallelOptions { MaxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads },
+                    new ParallelOptions
+                    {
+                        MaxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads
+                    },
                     (range, loopState) =>
                     {
                         for (var i = range.Item1; i < range.Item2; i++)
@@ -75,8 +110,15 @@ namespace MathNet.Numerics.Threading
                             body(i);
                         }
                     });
-            }
 #endif
+            }
+            else
+            {
+                for (var index = fromInclusive; index < toExclusive; index++)
+                {
+                    body(index);
+                }
+            }
         }
 
 /*        /// <summary>