diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj
index 12b2677e..b1df4e62 100644
--- a/src/FSharp/FSharp.fsproj
+++ b/src/FSharp/FSharp.fsproj
@@ -35,12 +35,12 @@
   </PropertyGroup>
   <ItemGroup>
     <Reference Include="FSharp.Core" />
-    <Reference Include="FSharp.Core, Version=2.0.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a" />
     <Reference Include="mscorlib" />
     <Reference Include="System" />
     <Reference Include="System.Core">
       <RequiredTargetFramework>3.5</RequiredTargetFramework>
     </Reference>
+    <Reference Include="System.Numerics" />
   </ItemGroup>
   <ItemGroup>
     <ProjectReference Include="..\Numerics\Numerics.csproj">
diff --git a/src/FSharpUnitTests/FSharpUnitTests.fsproj b/src/FSharpUnitTests/FSharpUnitTests.fsproj
index 54ac3ec6..2080d67e 100644
--- a/src/FSharpUnitTests/FSharpUnitTests.fsproj
+++ b/src/FSharpUnitTests/FSharpUnitTests.fsproj
@@ -30,6 +30,15 @@
     <DefineConstants>TRACE</DefineConstants>
     <WarningLevel>3</WarningLevel>
   </PropertyGroup>
+  <Import Project="$(MSBuildExtensionsPath32)\FSharp\1.0\Microsoft.FSharp.Targets" Condition="!Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <Import Project="$(MSBuildExtensionsPath32)\..\Microsoft F#\v4.0\Microsoft.FSharp.Targets" Condition="Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <ItemGroup>
+    <Compile Include="FsUnit.fs" />
+    <Compile Include="Program.fs" />
+    <None Include="App.config">
+      <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+    </None>
+  </ItemGroup>
   <ItemGroup>
     <Reference Include="FSharp.Core" />
     <Reference Include="mscorlib" />
@@ -37,12 +46,6 @@
     <Reference Include="System.Core">
       <RequiredTargetFramework>3.5</RequiredTargetFramework>
     </Reference>
-  </ItemGroup>
-  <ItemGroup>
-    <Compile Include="FsUnit.fs" />
-    <Compile Include="Program.fs" />
-  </ItemGroup>
-  <ItemGroup>
     <ProjectReference Include="..\FSharp\FSharp.fsproj">
       <Name>FSharp</Name>
       <Project>{37e8e802-a354-4114-bfc1-6e1357da605b}</Project>
@@ -53,14 +56,8 @@
       <Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
       <Private>True</Private>
     </ProjectReference>
+    <Reference Include="System.Numerics" />
   </ItemGroup>
-  <ItemGroup>
-    <None Include="App.config">
-      <CopyToOutputDirectory>Always</CopyToOutputDirectory>
-    </None>
-  </ItemGroup>
-  <Import Project="$(MSBuildExtensionsPath32)\FSharp\1.0\Microsoft.FSharp.Targets" Condition="!Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
-  <Import Project="$(MSBuildExtensionsPath32)\..\Microsoft F#\v4.0\Microsoft.FSharp.Targets" Condition="Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
 	     Other similar extension points exist, see Microsoft.Common.targets.
 	<Target Name="BeforeBuild">
diff --git a/src/MathNet.Numerics.5.0.ReSharper b/src/MathNet.Numerics.5.0.ReSharper
index 58e0d838..ca6c2a14 100644
--- a/src/MathNet.Numerics.5.0.ReSharper
+++ b/src/MathNet.Numerics.5.0.ReSharper
@@ -32,9 +32,12 @@
           <Item>unsafe</Item>
           <Item>volatile</Item>
         </MODIFIERS_ORDER>
+        <PLACE_ABSTRACT_ACCESSORHOLDER_ON_SINGLE_LINE>False</PLACE_ABSTRACT_ACCESSORHOLDER_ON_SINGLE_LINE>
         <PLACE_FIELD_ATTRIBUTE_ON_SAME_LINE>False</PLACE_FIELD_ATTRIBUTE_ON_SAME_LINE>
+        <PLACE_SIMPLE_ACCESSOR_ON_SINGLE_LINE>False</PLACE_SIMPLE_ACCESSOR_ON_SINGLE_LINE>
         <PLACE_TYPE_CONSTRAINTS_ON_SAME_LINE>False</PLACE_TYPE_CONSTRAINTS_ON_SAME_LINE>
         <PLACE_WHILE_ON_NEW_LINE>True</PLACE_WHILE_ON_NEW_LINE>
+        <REDUNDANT_THIS_QUALIFIER_STYLE>ALWAYS_USE</REDUNDANT_THIS_QUALIFIER_STYLE>
         <SPACE_AFTER_TYPECAST_PARENTHESES>False</SPACE_AFTER_TYPECAST_PARENTHESES>
         <SPACE_AROUND_MULTIPLICATIVE_OP>True</SPACE_AROUND_MULTIPLICATIVE_OP>
         <SPACE_BEFORE_SIZEOF_PARENTHESES>False</SPACE_BEFORE_SIZEOF_PARENTHESES>
@@ -45,6 +48,7 @@
         <WRAP_ARGUMENTS_STYLE>CHOP_IF_LONG</WRAP_ARGUMENTS_STYLE>
         <WRAP_BEFORE_BINARY_OPSIGN>True</WRAP_BEFORE_BINARY_OPSIGN>
         <WRAP_BEFORE_FIRST_TYPE_PARAMETER_CONSTRAINT>True</WRAP_BEFORE_FIRST_TYPE_PARAMETER_CONSTRAINT>
+        <WRAP_LINES>False</WRAP_LINES>
         <WRAP_OBJECT_AND_COLLECTION_INITIALIZER_STYLE>CHOP_ALWAYS</WRAP_OBJECT_AND_COLLECTION_INITIALIZER_STYLE>
         <WRAP_PARAMETERS_STYLE>CHOP_IF_LONG</WRAP_PARAMETERS_STYLE>
       </FormatSettings>
@@ -68,6 +72,10 @@
         <PredefinedRule Inspect="True" Prefix="_" Suffix="" Style="aaBb" ElementKind="StaticReadonly" />
         <PredefinedRule Inspect="True" Prefix="" Suffix="" Style="AaBb" ElementKind="EnumMember" />
         <PredefinedRule Inspect="True" Prefix="" Suffix="" Style="AaBb" ElementKind="Other" />
+        <PredefinedRule Inspect="True" Prefix="_" Suffix="" Style="aaBb" ElementKind="PrivateInstanceFields" />
+        <PredefinedRule Inspect="True" Prefix="_" Suffix="" Style="aaBb" ElementKind="PrivateStaticFields" />
+        <PredefinedRule Inspect="True" Prefix="" Suffix="" Style="AaBb" ElementKind="PrivateConstants" />
+        <PredefinedRule Inspect="True" Prefix="" Suffix="" Style="AaBb" ElementKind="PrivateStaticReadonly" />
       </Naming2>
     </CSharp>
     <VB>
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
index 2b1db7ea..f29344f3 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
@@ -24,9 +24,10 @@
 namespace MathNet.Numerics.Algorithms.LinearAlgebra
 {
     using System;
+    using System.Collections.Concurrent;
     using System.Numerics;
+    using System.Threading.Tasks;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// The managed linear algebra provider.
@@ -66,11 +67,27 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (alpha == 1.0)
             {
-                Parallel.For(0, y.Length, i => y[i] += x[i]);
+                Parallel.ForEach(
+                    Partitioner.Create(0, y.Length), 
+                    (range, loopState) =>
+                    {
+                        for (var i = range.Item1; i < range.Item2; i++)
+                        {
+                            y[i] += x[i];
+                        }
+                    });
             }
             else
             {
-                Parallel.For(0, y.Length, i => y[i] += alpha * x[i]);
+                Parallel.ForEach(
+                    Partitioner.Create(0, y.Length), 
+                    (range, loopState) =>
+                    {
+                        for (var i = range.Item1; i < range.Item2; i++)
+                        {
+                            y[i] += alpha * x[i];
+                        }
+                    });
             }
         }
 
@@ -92,7 +109,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 return;
             }
 
-            Parallel.For(0, x.Length, i => x[i] = alpha * x[i]);
+            Parallel.ForEach(
+                Partitioner.Create(0, x.Length), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        x[i] = alpha * x[i];
+                    }
+                });
         }
 
         /// <summary>
@@ -121,19 +146,26 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             var sync = new object();
             var d = 0.0;
-            Parallel.For(
-                0,
-                y.Length,
-                () => 0.0,
-                (i, localData) => localData += y[i] * x[i],
+
+            Parallel.ForEach(
+                Partitioner.Create(0, y.Length), 
+                () => 0.0, 
+                (range, loopState, localData) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        localData += y[i] * x[i];
+                    }
+
+                    return localData;
+                }, 
                 localResult =>
                 {
                     lock (sync)
                     {
                         d += localResult;
                     }
-                }
-                );
+                });
 
             return d;
         }
@@ -169,8 +201,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
-            
-            Parallel.For(0, y.Length, i => result[i] = x[i] + y[i]);
+
+            Parallel.ForEach(
+                Partitioner.Create(0, y.Length), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        result[i] = x[i] + y[i];
+                    }
+                });
         }
 
         /// <summary>
@@ -205,7 +245,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            Parallel.For(0, y.Length, i => result[i] = x[i] - y[i]);
+            Parallel.ForEach(
+                Partitioner.Create(0, y.Length), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        result[i] = x[i] - y[i];
+                    }
+                });
         }
 
         /// <summary>
@@ -239,8 +287,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
-
-            Parallel.For(0, y.Length, i => result[i] = x[i] * y[i]);
+            
+            Parallel.ForEach(
+                Partitioner.Create(0, y.Length), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        result[i] = x[i] * y[i];
+                    }
+                });
         }
 
         public double MatrixNorm(Norm norm, double[] matrix)
@@ -272,6 +328,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 throw new ArgumentNullException("x");
             }
+
             if (y == null)
             {
                 throw new ArgumentNullException("y");
@@ -308,7 +365,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             double[] xdata;
             if (ReferenceEquals(x, result))
             {
-                xdata = (double[]) x.Clone();
+                xdata = (double[])x.Clone();
             }
             else
             {
@@ -318,7 +375,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             double[] ydata;
             if (ReferenceEquals(y, result))
             {
-                ydata = (double[]) y.Clone();
+                ydata = (double[])y.Clone();
             }
             else
             {
@@ -329,7 +386,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             // TODO - For small matrices we should get rid of the parallelism because of startup costs.
             // Perhaps the following implementations would be a good one
             // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
-            MatrixMultiplyWithUpdate(Transpose.DontTranspose,Transpose.DontTranspose,1.0,x,xRows,xColumns,y,yRows,yColumns,0.0,result);
+            this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, xdata, xRows, xColumns, ydata, yRows, yColumns, 0.0, result);
         }
 
         /// <summary>
@@ -347,11 +404,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
         /// <param name="c">The c matrix.</param>
         public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, 
-            int aRows, int aColumns, double[] b, int bRows, int bColumns, double beta, double[] c)
+                                             int aRows, int aColumns, double[] b, int bRows, int bColumns, double beta, double[] c)
         {
             // Choose nonsensical values for the number of rows in c; fill them in depending
             // on the operations on a and b.
-            int cRows = -1;
+            var cRows = -1;
 
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (a == null)
@@ -427,7 +484,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 return;
             }
 
-
             // Check whether we will be overwriting any of our inputs and make copies if necessary.
             // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
             // as result, we can do it on a row wise basis. We should investigate this.
@@ -457,142 +513,190 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 {
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
-                        Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, aColumns), 
+                            (range, loopState) =>
                             {
-                                int iIndex = i * aRows;
-                                double s = 0;
-                                for (int l = 0; l != bColumns; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
+                                    var jIndex = j * cRows;
+                                    for (var i = 0; i != bRows; i++)
+                                    {
+                                        var iIndex = i * aRows;
+                                        double s = 0;
+                                        for (var l = 0; l != bColumns; l++)
+                                        {
+                                            s += adata[iIndex + l] * bdata[l * bRows + j];
+                                        }
+
+                                        c[jIndex + i] = s;
+                                    }
                                 }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                            });
                     }
                     else if ((int)transposeA > 111)
                     {
-                        Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, bColumns), 
+                            (range, loopState) =>
                             {
-                                int iIndex = i * aRows;
-                                double s = 0;
-                                for (int l = 0; l != aRows; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
+                                    var jcIndex = j * cRows;
+                                    var jbIndex = j * bRows;
+                                    for (var i = 0; i != aColumns; i++)
+                                    {
+                                        var iIndex = i * aRows;
+                                        double s = 0;
+                                        for (var l = 0; l != aRows; l++)
+                                        {
+                                            s += adata[iIndex + l] * bdata[jbIndex + l];
+                                        }
+
+                                        c[jcIndex + i] = s;
+                                    }
                                 }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                            });
                     }
                     else if ((int)transposeB > 111)
                     {
-                        Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, bRows), 
+                            (range, loopState) =>
                             {
-                                double s = 0;
-                                for (int l = 0; l != aColumns; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                    var jIndex = j * cRows;
+                                    for (var i = 0; i != aRows; i++)
+                                    {
+                                        double s = 0;
+                                        for (var l = 0; l != aColumns; l++)
+                                        {
+                                            s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                        }
+
+                                        c[jIndex + i] = s;
+                                    }
                                 }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                            });
                     }
                     else
                     {
-                        Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, bColumns), 
+                            (range, loopState) =>
                             {
-                                double s = 0;
-                                for (int l = 0; l != aColumns; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                    var jcIndex = j * cRows;
+                                    var jbIndex = j * bRows;
+                                    for (var i = 0; i != aRows; i++)
+                                    {
+                                        double s = 0;
+                                        for (var l = 0; l != aColumns; l++)
+                                        {
+                                            s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                        }
+
+                                        c[jcIndex + i] = s;
+                                    }
                                 }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                            });
                     }
                 }
                 else
                 {
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
-                        Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, aColumns), 
+                            (range, loopState) =>
                             {
-                                int iIndex = i * aRows;
-                                double s = 0;
-                                for (int l = 0; l != bColumns; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
+                                    var jIndex = j * cRows;
+                                    for (var i = 0; i != bRows; i++)
+                                    {
+                                        var iIndex = i * aRows;
+                                        double s = 0;
+                                        for (var l = 0; l != bColumns; l++)
+                                        {
+                                            s += adata[iIndex + l] * bdata[l * bRows + j];
+                                        }
+
+                                        c[jIndex + i] = c[jIndex + i] * beta + s;
+                                    }
                                 }
-                                c[jIndex + i] = c[jIndex + i] * beta + s;
-                            }
-                        });
+                            });
                     }
                     else if ((int)transposeA > 111)
                     {
-                        Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, bColumns), 
+                            (range, loopState) =>
                             {
-                                int iIndex = i * aRows;
-                                double s = 0;
-                                for (int l = 0; l != aRows; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
+                                    var jcIndex = j * cRows;
+                                    var jbIndex = j * bRows;
+                                    for (var i = 0; i != aColumns; i++)
+                                    {
+                                        var iIndex = i * aRows;
+                                        double s = 0;
+                                        for (var l = 0; l != aRows; l++)
+                                        {
+                                            s += adata[iIndex + l] * bdata[jbIndex + l];
+                                        }
+
+                                        c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                    }
                                 }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                            });
                     }
                     else if ((int)transposeB > 111)
                     {
-                        Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, bRows), 
+                            (range, loopState) =>
                             {
-                                double s = 0;
-                                for (int l = 0; l != aColumns; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                    var jIndex = j * cRows;
+                                    for (var i = 0; i != aRows; i++)
+                                    {
+                                        double s = 0;
+                                        for (var l = 0; l != aColumns; l++)
+                                        {
+                                            s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                        }
+
+                                        c[jIndex + i] = s + c[jIndex + i] * beta;
+                                    }
                                 }
-                                c[jIndex + i] = s + c[jIndex + i] * beta;
-                            }
-                        });
+                            });
                     }
                     else
                     {
-                        Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
+                        Parallel.ForEach(
+                            Partitioner.Create(0, bColumns), 
+                            (range, loopState) =>
                             {
-                                double s = 0;
-                                for (int l = 0; l != aColumns; l++)
+                                for (var j = range.Item1; j < range.Item2; j++)
                                 {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                    var jcIndex = j * cRows;
+                                    var jbIndex = j * bRows;
+                                    for (var i = 0; i != aRows; i++)
+                                    {
+                                        double s = 0;
+                                        for (var l = 0; l != aColumns; l++)
+                                        {
+                                            s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                        }
+
+                                        c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                    }
                                 }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                            });
                     }
                 }
             }
@@ -600,71 +704,95 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 if ((int)transposeA > 111 && (int)transposeB > 111)
                 {
-                    Parallel.For(0, aColumns, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != bRows; i++)
+                    Parallel.ForEach(
+                        Partitioner.Create(0, aColumns), 
+                        (range, loopState) =>
                         {
-                            int iIndex = i * aRows;
-                            double s = 0;
-                            for (int l = 0; l != bColumns; l++)
+                            for (var j = range.Item1; j < range.Item2; j++)
                             {
-                                s += adata[iIndex + l] * bdata[l * bRows + j];
+                                var jIndex = j * cRows;
+                                for (var i = 0; i != bRows; i++)
+                                {
+                                    var iIndex = i * aRows;
+                                    double s = 0;
+                                    for (var l = 0; l != bColumns; l++)
+                                    {
+                                        s += adata[iIndex + l] * bdata[l * bRows + j];
+                                    }
+
+                                    c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
+                                }
                             }
-                            c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
-                        }
-                    });
+                        });
                 }
                 else if ((int)transposeA > 111)
                 {
-                    Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aColumns; i++)
+                    Parallel.ForEach(
+                        Partitioner.Create(0, bColumns), 
+                        (range, loopState) =>
                         {
-                            int iIndex = i * aRows;
-                            double s = 0;
-                            for (int l = 0; l != aRows; l++)
+                            for (var j = range.Item1; j < range.Item2; j++)
                             {
-                                s += adata[iIndex + l] * bdata[jbIndex + l];
+                                var jcIndex = j * cRows;
+                                var jbIndex = j * bRows;
+                                for (var i = 0; i != aColumns; i++)
+                                {
+                                    var iIndex = i * aRows;
+                                    double s = 0;
+                                    for (var l = 0; l != aRows; l++)
+                                    {
+                                        s += adata[iIndex + l] * bdata[jbIndex + l];
+                                    }
+
+                                    c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                }
                             }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                        });
                 }
                 else if ((int)transposeB > 111)
                 {
-                    Parallel.For(0, bRows, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != aRows; i++)
+                    Parallel.ForEach(
+                        Partitioner.Create(0, bRows), 
+                        (range, loopState) =>
                         {
-                            double s = 0;
-                            for (int l = 0; l != aColumns; l++)
+                            for (var j = range.Item1; j < range.Item2; j++)
                             {
-                                s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                var jIndex = j * cRows;
+                                for (var i = 0; i != aRows; i++)
+                                {
+                                    double s = 0;
+                                    for (var l = 0; l != aColumns; l++)
+                                    {
+                                        s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                    }
+
+                                    c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
+                                }
                             }
-                            c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
-                        }
-                    });
+                        });
                 }
                 else
                 {
-                    Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aRows; i++)
+                    Parallel.ForEach(
+                        Partitioner.Create(0, bColumns), 
+                        (range, loopState) =>
                         {
-                            double s = 0;
-                            for (int l = 0; l != aColumns; l++)
+                            for (var j = range.Item1; j < range.Item2; j++)
                             {
-                                s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                var jcIndex = j * cRows;
+                                var jbIndex = j * bRows;
+                                for (var i = 0; i != aRows; i++)
+                                {
+                                    double s = 0;
+                                    for (var l = 0; l != aColumns; l++)
+                                    {
+                                        s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                    }
+
+                                    c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                }
                             }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                        });
                 }
             }
         }
@@ -735,6 +863,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     {
                         s += a[i * order + k] * a[i * order + j];
                     }
+
                     var tmp = k * order;
                     index = tmp + j;
                     a[index] = s = (a[index] - s) / a[tmp + k];
@@ -904,9 +1033,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            float d = 0.0F;
+            var d = 0.0F;
 
-            for (int i = 0; i < y.Length; i++)
+            for (var i = 0; i < y.Length; i++)
             {
                 d += y[i] * x[i];
             }
@@ -1048,6 +1177,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 throw new ArgumentNullException("x");
             }
+
             if (y == null)
             {
                 throw new ArgumentNullException("y");
@@ -1084,7 +1214,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             float[] xdata;
             if (ReferenceEquals(x, result))
             {
-                xdata = (float[]) x.Clone();
+                xdata = (float[])x.Clone();
             }
             else
             {
@@ -1094,7 +1224,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             float[] ydata;
             if (ReferenceEquals(y, result))
             {
-                ydata = (float[]) y.Clone();
+                ydata = (float[])y.Clone();
             }
             else
             {
@@ -1105,8 +1235,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             // TODO - For small matrices we should get rid of the parallelism because of startup costs.
             // Perhaps the following implementations would be a good one
             // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
-            MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, xRows, xColumns, y, yRows, yColumns, 0.0f, result);
-
+            this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, xRows, xColumns, y, yRows, yColumns, 0.0f, result);
         }
 
         /// <summary>
@@ -1124,12 +1253,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
         /// <param name="c">The c matrix.</param>
         public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, 
-            int aRows, int aColumns, float[] b, int bRows, int bColumns, float beta, float[] c)
+                                             int aRows, int aColumns, float[] b, int bRows, int bColumns, float beta, float[] c)
         {
             // Choose nonsensical values for the number of rows and columns in c; fill them in depending
             // on the operations on a and b.
-            int cRows = -1;
-            int cColumns = -1;
+            var cRows = -1;
+            var cColumns = -1;
 
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (a == null)
@@ -1240,70 +1369,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
                         Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
-                            {
-                                int iIndex = i * aRows;
-                                float s = 0;
-                                for (int l = 0; l != bColumns; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jIndex = j * cRows;
+                                                      for (var i = 0; i != bRows; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          float s = 0;
+                                                          for (var l = 0; l != bColumns; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                          }
+
+                                                          c[jIndex + i] = s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeA > 111)
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
-                            {
-                                int iIndex = i * aRows;
-                                float s = 0;
-                                for (int l = 0; l != aRows; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aColumns; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          float s = 0;
+                                                          for (var l = 0; l != aRows; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeB > 111)
                     {
                         Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                float s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                                               {
+                                                   var jIndex = j * cRows;
+                                                   for (var i = 0; i != aRows; i++)
+                                                   {
+                                                       float s = 0;
+                                                       for (var l = 0; l != aColumns; l++)
+                                                       {
+                                                           s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                       }
+
+                                                       c[jIndex + i] = s;
+                                                   }
+                                               });
                     }
                     else
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                float s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aRows; i++)
+                                                      {
+                                                          float s = 0;
+                                                          for (var l = 0; l != aColumns; l++)
+                                                          {
+                                                              s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s;
+                                                      }
+                                                  });
                     }
                 }
                 else
@@ -1311,70 +1444,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
                         Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
-                            {
-                                int iIndex = i * aRows;
-                                float s = 0;
-                                for (int l = 0; l != bColumns; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = c[jIndex + i] * beta + s;
-                            }
-                        });
+                                                  {
+                                                      var jIndex = j * cRows;
+                                                      for (var i = 0; i != bRows; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          float s = 0;
+                                                          for (var l = 0; l != bColumns; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                          }
+
+                                                          c[jIndex + i] = c[jIndex + i] * beta + s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeA > 111)
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
-                            {
-                                int iIndex = i * aRows;
-                                float s = 0;
-                                for (int l = 0; l != aRows; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aColumns; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          float s = 0;
+                                                          for (var l = 0; l != aRows; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeB > 111)
                     {
                         Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                float s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s + c[jIndex + i] * beta;
-                            }
-                        });
+                                               {
+                                                   var jIndex = j * cRows;
+                                                   for (var i = 0; i != aRows; i++)
+                                                   {
+                                                       float s = 0;
+                                                       for (var l = 0; l != aColumns; l++)
+                                                       {
+                                                           s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                       }
+
+                                                       c[jIndex + i] = s + c[jIndex + i] * beta;
+                                                   }
+                                               });
                     }
                     else
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                float s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aRows; i++)
+                                                      {
+                                                          float s = 0;
+                                                          for (var l = 0; l != aColumns; l++)
+                                                          {
+                                                              s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                                      }
+                                                  });
                     }
                 }
             }
@@ -1383,70 +1520,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 if ((int)transposeA > 111 && (int)transposeB > 111)
                 {
                     Parallel.For(0, aColumns, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != bRows; i++)
-                        {
-                            int iIndex = i * aRows;
-                            float s = 0;
-                            for (int l = 0; l != bColumns; l++)
-                            {
-                                s += adata[iIndex + l] * bdata[l * bRows + j];
-                            }
-                            c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
-                        }
-                    });
+                                              {
+                                                  var jIndex = j * cRows;
+                                                  for (var i = 0; i != bRows; i++)
+                                                  {
+                                                      var iIndex = i * aRows;
+                                                      float s = 0;
+                                                      for (var l = 0; l != bColumns; l++)
+                                                      {
+                                                          s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                      }
+
+                                                      c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
+                                                  }
+                                              });
                 }
                 else if ((int)transposeA > 111)
                 {
                     Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aColumns; i++)
-                        {
-                            int iIndex = i * aRows;
-                            float s = 0;
-                            for (int l = 0; l != aRows; l++)
-                            {
-                                s += adata[iIndex + l] * bdata[jbIndex + l];
-                            }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                                              {
+                                                  var jcIndex = j * cRows;
+                                                  var jbIndex = j * bRows;
+                                                  for (var i = 0; i != aColumns; i++)
+                                                  {
+                                                      var iIndex = i * aRows;
+                                                      float s = 0;
+                                                      for (var l = 0; l != aRows; l++)
+                                                      {
+                                                          s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                      }
+
+                                                      c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                                  }
+                                              });
                 }
                 else if ((int)transposeB > 111)
                 {
                     Parallel.For(0, bRows, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != aRows; i++)
-                        {
-                            float s = 0;
-                            for (int l = 0; l != aColumns; l++)
-                            {
-                                s += adata[l * aRows + i] * bdata[l * bRows + j];
-                            }
-                            c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
-                        }
-                    });
+                                           {
+                                               var jIndex = j * cRows;
+                                               for (var i = 0; i != aRows; i++)
+                                               {
+                                                   float s = 0;
+                                                   for (var l = 0; l != aColumns; l++)
+                                                   {
+                                                       s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                   }
+
+                                                   c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
+                                               }
+                                           });
                 }
                 else
                 {
                     Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aRows; i++)
-                        {
-                            float s = 0;
-                            for (int l = 0; l != aColumns; l++)
-                            {
-                                s += adata[l * aRows + i] * bdata[jbIndex + l];
-                            }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                                              {
+                                                  var jcIndex = j * cRows;
+                                                  var jbIndex = j * bRows;
+                                                  for (var i = 0; i != aRows; i++)
+                                                  {
+                                                      float s = 0;
+                                                      for (var l = 0; l != aColumns; l++)
+                                                      {
+                                                          s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                      }
+
+                                                      c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                                  }
+                                              });
                 }
             }
         }
@@ -1505,33 +1646,36 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
         public void CholeskyFactor(float[] a, int order)
         {
-            float[] factor = new float[a.Length];
+            var factor = new float[a.Length];
 
-            for (int j = 0; j < order; j++)
+            for (var j = 0; j < order; j++)
             {
-                float d = 0.0F;
+                var d = 0.0F;
                 int index;
-                for (int k = 0; k < j; k++)
+                for (var k = 0; k < j; k++)
                 {
-                    float s = 0.0F;
+                    var s = 0.0F;
                     int i;
                     for (i = 0; i < k; i++)
                     {
                         s += factor[i * order + k] * factor[i * order + j];
                     }
-                    int tmp = k * order;
+
+                    var tmp = k * order;
                     index = tmp + j;
                     factor[index] = s = (a[index] - s) / factor[tmp + k];
                     d += s * s;
                 }
+
                 index = j * order + j;
                 d = a[index] - d;
                 if (d <= 0.0F)
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
-                factor[index] = (float) System.Math.Sqrt(d);
-                for (int k = j + 1; k < order; k++)
+
+                factor[index] = (float)Math.Sqrt(d);
+                for (var k = j + 1; k < order; k++)
                 {
                     factor[k * order + j] = 0.0F;
                 }
@@ -1688,9 +1832,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            Complex d = new Complex(0.0, 0.0);
+            var d = new Complex(0.0, 0.0);
 
-            for (int i = 0; i < y.Length; i++)
+            for (var i = 0; i < y.Length; i++)
             {
                 d += y[i] * x[i];
             }
@@ -1812,7 +1956,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         {
             throw new NotImplementedException();
         }
-        
+
         /// <summary>
         /// Multiples two matrices. <c>result = x * y</c>
         /// </summary>
@@ -1832,6 +1976,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 throw new ArgumentNullException("x");
             }
+
             if (y == null)
             {
                 throw new ArgumentNullException("y");
@@ -1868,7 +2013,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             Complex[] xdata;
             if (ReferenceEquals(x, result))
             {
-                xdata = (Complex[]) x.Clone();
+                xdata = (Complex[])x.Clone();
             }
             else
             {
@@ -1878,7 +2023,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             Complex[] ydata;
             if (ReferenceEquals(y, result))
             {
-                ydata = (Complex[]) y.Clone();
+                ydata = (Complex[])y.Clone();
             }
             else
             {
@@ -1889,7 +2034,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             // TODO - For small matrices we should get rid of the parallelism because of startup costs.
             // Perhaps the following implementations would be a good one
             // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
-            MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, xRows, xColumns, y, yRows, yColumns, Complex.Zero, result);
+            this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, xRows, xColumns, y, yRows, yColumns, Complex.Zero, result);
         }
 
         /// <summary>
@@ -1907,12 +2052,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
         /// <param name="c">The c matrix.</param>
         public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, 
-            int aRows, int aColumns, Complex[] b, int bRows, int bColumns, Complex beta, Complex[] c)
+                                             int aRows, int aColumns, Complex[] b, int bRows, int bColumns, Complex beta, Complex[] c)
         {
             // Choose nonsensical values for the number of rows and columns in c; fill them in depending
             // on the operations on a and b.
-            int cRows = -1;
-            int cColumns = -1;
+            var cRows = -1;
+            var cColumns = -1;
 
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (a == null)
@@ -2023,70 +2168,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
                         Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex s = 0;
-                                for (int l = 0; l != bColumns; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jIndex = j * cRows;
+                                                      for (var i = 0; i != bRows; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex s = 0;
+                                                          for (var l = 0; l != bColumns; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                          }
+
+                                                          c[jIndex + i] = s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeA > 111)
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex s = 0;
-                                for (int l = 0; l != aRows; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aColumns; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex s = 0;
+                                                          for (var l = 0; l != aRows; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeB > 111)
                     {
                         Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                                               {
+                                                   var jIndex = j * cRows;
+                                                   for (var i = 0; i != aRows; i++)
+                                                   {
+                                                       Complex s = 0;
+                                                       for (var l = 0; l != aColumns; l++)
+                                                       {
+                                                           s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                       }
+
+                                                       c[jIndex + i] = s;
+                                                   }
+                                               });
                     }
                     else
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aRows; i++)
+                                                      {
+                                                          Complex s = 0;
+                                                          for (var l = 0; l != aColumns; l++)
+                                                          {
+                                                              s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s;
+                                                      }
+                                                  });
                     }
                 }
                 else
@@ -2094,70 +2243,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
                         Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex s = 0;
-                                for (int l = 0; l != bColumns; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = c[jIndex + i] * beta + s;
-                            }
-                        });
+                                                  {
+                                                      var jIndex = j * cRows;
+                                                      for (var i = 0; i != bRows; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex s = 0;
+                                                          for (var l = 0; l != bColumns; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                          }
+
+                                                          c[jIndex + i] = c[jIndex + i] * beta + s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeA > 111)
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex s = 0;
-                                for (int l = 0; l != aRows; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aColumns; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex s = 0;
+                                                          for (var l = 0; l != aRows; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeB > 111)
                     {
                         Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s + c[jIndex + i] * beta;
-                            }
-                        });
+                                               {
+                                                   var jIndex = j * cRows;
+                                                   for (var i = 0; i != aRows; i++)
+                                                   {
+                                                       Complex s = 0;
+                                                       for (var l = 0; l != aColumns; l++)
+                                                       {
+                                                           s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                       }
+
+                                                       c[jIndex + i] = s + c[jIndex + i] * beta;
+                                                   }
+                                               });
                     }
                     else
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aRows; i++)
+                                                      {
+                                                          Complex s = 0;
+                                                          for (var l = 0; l != aColumns; l++)
+                                                          {
+                                                              s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                                      }
+                                                  });
                     }
                 }
             }
@@ -2166,70 +2319,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 if ((int)transposeA > 111 && (int)transposeB > 111)
                 {
                     Parallel.For(0, aColumns, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != bRows; i++)
-                        {
-                            int iIndex = i * aRows;
-                            Complex s = 0;
-                            for (int l = 0; l != bColumns; l++)
-                            {
-                                s += adata[iIndex + l] * bdata[l * bRows + j];
-                            }
-                            c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
-                        }
-                    });
+                                              {
+                                                  var jIndex = j * cRows;
+                                                  for (var i = 0; i != bRows; i++)
+                                                  {
+                                                      var iIndex = i * aRows;
+                                                      Complex s = 0;
+                                                      for (var l = 0; l != bColumns; l++)
+                                                      {
+                                                          s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                      }
+
+                                                      c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
+                                                  }
+                                              });
                 }
                 else if ((int)transposeA > 111)
                 {
                     Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aColumns; i++)
-                        {
-                            int iIndex = i * aRows;
-                            Complex s = 0;
-                            for (int l = 0; l != aRows; l++)
-                            {
-                                s += adata[iIndex + l] * bdata[jbIndex + l];
-                            }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                                              {
+                                                  var jcIndex = j * cRows;
+                                                  var jbIndex = j * bRows;
+                                                  for (var i = 0; i != aColumns; i++)
+                                                  {
+                                                      var iIndex = i * aRows;
+                                                      Complex s = 0;
+                                                      for (var l = 0; l != aRows; l++)
+                                                      {
+                                                          s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                      }
+
+                                                      c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                                  }
+                                              });
                 }
                 else if ((int)transposeB > 111)
                 {
                     Parallel.For(0, bRows, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != aRows; i++)
-                        {
-                            Complex s = 0;
-                            for (int l = 0; l != aColumns; l++)
-                            {
-                                s += adata[l * aRows + i] * bdata[l * bRows + j];
-                            }
-                            c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
-                        }
-                    });
+                                           {
+                                               var jIndex = j * cRows;
+                                               for (var i = 0; i != aRows; i++)
+                                               {
+                                                   Complex s = 0;
+                                                   for (var l = 0; l != aColumns; l++)
+                                                   {
+                                                       s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                   }
+
+                                                   c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
+                                               }
+                                           });
                 }
                 else
                 {
                     Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aRows; i++)
-                        {
-                            Complex s = 0;
-                            for (int l = 0; l != aColumns; l++)
-                            {
-                                s += adata[l * aRows + i] * bdata[jbIndex + l];
-                            }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                                              {
+                                                  var jcIndex = j * cRows;
+                                                  var jbIndex = j * bRows;
+                                                  for (var i = 0; i != aRows; i++)
+                                                  {
+                                                      Complex s = 0;
+                                                      for (var l = 0; l != aColumns; l++)
+                                                      {
+                                                          s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                      }
+
+                                                      c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                                  }
+                                              });
                 }
             }
         }
@@ -2439,9 +2596,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            Complex32 d = new Complex32(0.0F, 0.0F);
+            var d = new Complex32(0.0F, 0.0F);
 
-            for (int i = 0; i < y.Length; i++)
+            for (var i = 0; i < y.Length; i++)
             {
                 d += y[i] * x[i];
             }
@@ -2583,6 +2740,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 throw new ArgumentNullException("x");
             }
+
             if (y == null)
             {
                 throw new ArgumentNullException("y");
@@ -2619,7 +2777,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             Complex32[] xdata;
             if (ReferenceEquals(x, result))
             {
-                xdata = (Complex32[]) x.Clone();
+                xdata = (Complex32[])x.Clone();
             }
             else
             {
@@ -2629,7 +2787,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             Complex32[] ydata;
             if (ReferenceEquals(y, result))
             {
-                ydata = (Complex32[]) y.Clone();
+                ydata = (Complex32[])y.Clone();
             }
             else
             {
@@ -2640,7 +2798,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             // TODO - For small matrices we should get rid of the parallelism because of startup costs.
             // Perhaps the following implementations would be a good one
             // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
-            MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, xRows, xColumns, y, yRows, yColumns, Complex32.Zero, result);
+            this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, xRows, xColumns, y, yRows, yColumns, Complex32.Zero, result);
         }
 
         /// <summary>
@@ -2658,12 +2816,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
         /// <param name="c">The c matrix.</param>
         public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, 
-            int aRows, int aColumns, Complex32[] b, int bRows, int bColumns, Complex32 beta, Complex32[] c)
+                                             int aRows, int aColumns, Complex32[] b, int bRows, int bColumns, Complex32 beta, Complex32[] c)
         {
             // Choose nonsensical values for the number of rows and columns in c; fill them in depending
             // on the operations on a and b.
-            int cRows = -1;
-            int cColumns = -1;
+            var cRows = -1;
+            var cColumns = -1;
 
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (a == null)
@@ -2774,70 +2932,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
                         Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex32 s = 0;
-                                for (int l = 0; l != bColumns; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jIndex = j * cRows;
+                                                      for (var i = 0; i != bRows; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex32 s = 0;
+                                                          for (var l = 0; l != bColumns; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                          }
+
+                                                          c[jIndex + i] = s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeA > 111)
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex32 s = 0;
-                                for (int l = 0; l != aRows; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aColumns; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex32 s = 0;
+                                                          for (var l = 0; l != aRows; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeB > 111)
                     {
                         Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex32 s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s;
-                            }
-                        });
+                                               {
+                                                   var jIndex = j * cRows;
+                                                   for (var i = 0; i != aRows; i++)
+                                                   {
+                                                       Complex32 s = 0;
+                                                       for (var l = 0; l != aColumns; l++)
+                                                       {
+                                                           s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                       }
+
+                                                       c[jIndex + i] = s;
+                                                   }
+                                               });
                     }
                     else
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex32 s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aRows; i++)
+                                                      {
+                                                          Complex32 s = 0;
+                                                          for (var l = 0; l != aColumns; l++)
+                                                          {
+                                                              s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s;
+                                                      }
+                                                  });
                     }
                 }
                 else
@@ -2845,70 +3007,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     if ((int)transposeA > 111 && (int)transposeB > 111)
                     {
                         Parallel.For(0, aColumns, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != bRows; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex32 s = 0;
-                                for (int l = 0; l != bColumns; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = c[jIndex + i] * beta + s;
-                            }
-                        });
+                                                  {
+                                                      var jIndex = j * cRows;
+                                                      for (var i = 0; i != bRows; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex32 s = 0;
+                                                          for (var l = 0; l != bColumns; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                          }
+
+                                                          c[jIndex + i] = c[jIndex + i] * beta + s;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeA > 111)
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aColumns; i++)
-                            {
-                                int iIndex = i * aRows;
-                                Complex32 s = 0;
-                                for (int l = 0; l != aRows; l++)
-                                {
-                                    s += adata[iIndex + l] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aColumns; i++)
+                                                      {
+                                                          var iIndex = i * aRows;
+                                                          Complex32 s = 0;
+                                                          for (var l = 0; l != aRows; l++)
+                                                          {
+                                                              s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                                      }
+                                                  });
                     }
                     else if ((int)transposeB > 111)
                     {
                         Parallel.For(0, bRows, j =>
-                        {
-                            int jIndex = j * cRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex32 s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[l * bRows + j];
-                                }
-                                c[jIndex + i] = s + c[jIndex + i] * beta;
-                            }
-                        });
+                                               {
+                                                   var jIndex = j * cRows;
+                                                   for (var i = 0; i != aRows; i++)
+                                                   {
+                                                       Complex32 s = 0;
+                                                       for (var l = 0; l != aColumns; l++)
+                                                       {
+                                                           s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                       }
+
+                                                       c[jIndex + i] = s + c[jIndex + i] * beta;
+                                                   }
+                                               });
                     }
                     else
                     {
                         Parallel.For(0, bColumns, j =>
-                        {
-                            int jcIndex = j * cRows;
-                            int jbIndex = j * bRows;
-                            for (int i = 0; i != aRows; i++)
-                            {
-                                Complex32 s = 0;
-                                for (int l = 0; l != aColumns; l++)
-                                {
-                                    s += adata[l * aRows + i] * bdata[jbIndex + l];
-                                }
-                                c[jcIndex + i] = s + c[jcIndex + i] * beta;
-                            }
-                        });
+                                                  {
+                                                      var jcIndex = j * cRows;
+                                                      var jbIndex = j * bRows;
+                                                      for (var i = 0; i != aRows; i++)
+                                                      {
+                                                          Complex32 s = 0;
+                                                          for (var l = 0; l != aColumns; l++)
+                                                          {
+                                                              s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                          }
+
+                                                          c[jcIndex + i] = s + c[jcIndex + i] * beta;
+                                                      }
+                                                  });
                     }
                 }
             }
@@ -2917,70 +3083,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 if ((int)transposeA > 111 && (int)transposeB > 111)
                 {
                     Parallel.For(0, aColumns, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != bRows; i++)
-                        {
-                            int iIndex = i * aRows;
-                            Complex32 s = 0;
-                            for (int l = 0; l != bColumns; l++)
-                            {
-                                s += adata[iIndex + l] * bdata[l * bRows + j];
-                            }
-                            c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
-                        }
-                    });
+                                              {
+                                                  var jIndex = j * cRows;
+                                                  for (var i = 0; i != bRows; i++)
+                                                  {
+                                                      var iIndex = i * aRows;
+                                                      Complex32 s = 0;
+                                                      for (var l = 0; l != bColumns; l++)
+                                                      {
+                                                          s += adata[iIndex + l] * bdata[l * bRows + j];
+                                                      }
+
+                                                      c[jIndex + i] = c[jIndex + i] * beta + alpha * s;
+                                                  }
+                                              });
                 }
                 else if ((int)transposeA > 111)
                 {
                     Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aColumns; i++)
-                        {
-                            int iIndex = i * aRows;
-                            Complex32 s = 0;
-                            for (int l = 0; l != aRows; l++)
-                            {
-                                s += adata[iIndex + l] * bdata[jbIndex + l];
-                            }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                                              {
+                                                  var jcIndex = j * cRows;
+                                                  var jbIndex = j * bRows;
+                                                  for (var i = 0; i != aColumns; i++)
+                                                  {
+                                                      var iIndex = i * aRows;
+                                                      Complex32 s = 0;
+                                                      for (var l = 0; l != aRows; l++)
+                                                      {
+                                                          s += adata[iIndex + l] * bdata[jbIndex + l];
+                                                      }
+
+                                                      c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                                  }
+                                              });
                 }
                 else if ((int)transposeB > 111)
                 {
                     Parallel.For(0, bRows, j =>
-                    {
-                        int jIndex = j * cRows;
-                        for (int i = 0; i != aRows; i++)
-                        {
-                            Complex32 s = 0;
-                            for (int l = 0; l != aColumns; l++)
-                            {
-                                s += adata[l * aRows + i] * bdata[l * bRows + j];
-                            }
-                            c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
-                        }
-                    });
+                                           {
+                                               var jIndex = j * cRows;
+                                               for (var i = 0; i != aRows; i++)
+                                               {
+                                                   Complex32 s = 0;
+                                                   for (var l = 0; l != aColumns; l++)
+                                                   {
+                                                       s += adata[l * aRows + i] * bdata[l * bRows + j];
+                                                   }
+
+                                                   c[jIndex + i] = alpha * s + c[jIndex + i] * beta;
+                                               }
+                                           });
                 }
                 else
                 {
                     Parallel.For(0, bColumns, j =>
-                    {
-                        int jcIndex = j * cRows;
-                        int jbIndex = j * bRows;
-                        for (int i = 0; i != aRows; i++)
-                        {
-                            Complex32 s = 0;
-                            for (int l = 0; l != aColumns; l++)
-                            {
-                                s += adata[l * aRows + i] * bdata[jbIndex + l];
-                            }
-                            c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
-                        }
-                    });
+                                              {
+                                                  var jcIndex = j * cRows;
+                                                  var jbIndex = j * bRows;
+                                                  for (var i = 0; i != aRows; i++)
+                                                  {
+                                                      Complex32 s = 0;
+                                                      for (var l = 0; l != aColumns; l++)
+                                                      {
+                                                          s += adata[l * aRows + i] * bdata[jbIndex + l];
+                                                      }
+
+                                                      c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta;
+                                                  }
+                                              });
                 }
             }
         }
diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs
index e762d4b3..15381230 100644
--- a/src/Numerics/Control.cs
+++ b/src/Numerics/Control.cs
@@ -29,7 +29,6 @@
 namespace MathNet.Numerics
 {
     using Algorithms.LinearAlgebra;
-    using Threading;
 
     /// <summary>
     /// Sets parameters for the library.
@@ -63,16 +62,6 @@ namespace MathNet.Numerics
         /// </value>
         public static bool ThreadSafeRandomNumberGenerators { get; set; }
 
-        /// <summary>
-        /// Gets or sets a value indicating how many parallel worker threads shall be used
-        /// when parallelization is applicable.
-        /// </summary>
-        public static int NumberOfParallelWorkerThreads
-        {
-            get { return ThreadQueue.ThreadCount; }
-            set { ThreadQueue.Start(value); }
-        }
-
         /// <summary>
         /// Gets or sets a value indicating whether parallelization shall be disabled globally.
         /// </summary>
diff --git a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Bluestein.cs b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Bluestein.cs
index cb70a23d..825ba0eb 100644
--- a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Bluestein.cs
+++ b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Bluestein.cs
@@ -30,9 +30,9 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
 {
     using System;
     using System.Numerics;
+    using System.Threading.Tasks;
     using NumberTheory;
-    using Threading;
-
+    
     /// <summary>
     /// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
     /// </summary>
@@ -71,7 +71,7 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
             Complex[] b = new Complex[m];
             Complex[] a = new Complex[m];
 
-            Parallel.Run(
+            Parallel.Invoke(
                 () =>
                 {
                     // Build and transform padded sequence b_k = exp(I*Pi*k^2/N)
diff --git a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs
index 1803aa37..cf9aadf9 100644
--- a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs
+++ b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs
@@ -1,9 +1,7 @@
 // <copyright file="DiscreteFourierTransform.Naive.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
-//
 // Copyright (c) 2009 Math.NET
-//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,10 +10,8 @@
 // copies of the Software, and to permit persons to whom the
 // Software is furnished to do so, subject to the following
 // conditions:
-//
 // The above copyright notice and this permission notice shall be
 // included in all copies or substantial portions of the Software.
-//
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
@@ -29,8 +25,9 @@
 namespace MathNet.Numerics.IntegralTransforms.Algorithms
 {
     using System;
+    using System.Collections.Concurrent;
     using System.Numerics;
-    using Threading;
+    using System.Threading.Tasks;
 
     /// <summary>
     /// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
@@ -45,23 +42,25 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
         /// <returns>Corresponding frequency-space vector.</returns>
         internal static Complex[] Naive(Complex[] samples, int exponentSign)
         {
-            double w0 = exponentSign * Constants.Pi2 / samples.Length;
+            var w0 = exponentSign * Constants.Pi2 / samples.Length;
             var spectrum = new Complex[samples.Length];
 
-            Parallel.For(
-                0,
-                samples.Length,
-                k =>
+            Parallel.ForEach(
+                Partitioner.Create(0, samples.Length), 
+                (range, loopState) =>
                 {
-                    double wk = w0 * k;
-                    Complex sum = Complex.Zero;
-                    for (int n = 0; n < samples.Length; n++)
+                    for (var k = range.Item1; k < range.Item2; k++)
                     {
-                        double w = n * wk;
-                        sum += samples[n] * new Complex(Math.Cos(w), Math.Sin(w));
-                    }
+                        var wk = w0 * k;
+                        var sum = Complex.Zero;
+                        for (var n = 0; n < samples.Length; n++)
+                        {
+                            var w = n * wk;
+                            sum += samples[n] * new Complex(Math.Cos(w), Math.Sin(w));
+                        }
 
-                    spectrum[k] = sum;
+                        spectrum[k] = sum;
+                    }
                 });
 
             return spectrum;
@@ -93,4 +92,4 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
             return timeSpace;
         }
     }
-}
+}
\ No newline at end of file
diff --git a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs
index 046bb40d..76107de1 100644
--- a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs
+++ b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs
@@ -1,9 +1,7 @@
 // <copyright file="DiscreteFourierTransform.RadixN.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
-//
 // Copyright (c) 2009 Math.NET
-//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,10 +10,8 @@
 // copies of the Software, and to permit persons to whom the
 // Software is furnished to do so, subject to the following
 // conditions:
-//
 // The above copyright notice and this permission notice shall be
 // included in all copies or substantial portions of the Software.
-//
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
@@ -29,10 +25,11 @@
 namespace MathNet.Numerics.IntegralTransforms.Algorithms
 {
     using System;
+    using System.Collections.Concurrent;
     using System.Numerics;
+    using System.Threading.Tasks;
     using NumberTheory;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
@@ -46,17 +43,17 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
         /// <param name="samples">Sample vector</param>
         private static void Radix2Reorder<T>(T[] samples)
         {
-            int j = 0;
-            for (int i = 0; i < samples.Length - 1; i++)
+            var j = 0;
+            for (var i = 0; i < samples.Length - 1; i++)
             {
                 if (i < j)
                 {
-                    T temp = samples[i];
+                    var temp = samples[i];
                     samples[i] = samples[j];
                     samples[j] = temp;
                 }
 
-                int m = samples.Length;
+                var m = samples.Length;
 
                 do
                 {
@@ -77,14 +74,14 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
         private static void Radix2Step(Complex[] samples, int exponentSign, int levelSize, int k)
         {
             // Twiddle Factor
-            double exponent = (exponentSign * k) * Constants.Pi / levelSize;
-            Complex w = new Complex(Math.Cos(exponent), Math.Sin(exponent));
+            var exponent = (exponentSign * k) * Constants.Pi / levelSize;
+            var w = new Complex(Math.Cos(exponent), Math.Sin(exponent));
 
-            int step = levelSize << 1;
-            for (int i = k; i < samples.Length; i += step)
+            var step = levelSize << 1;
+            for (var i = k; i < samples.Length; i += step)
             {
-                Complex ai = samples[i];
-                Complex t = w * samples[i + levelSize];
+                var ai = samples[i];
+                var t = w * samples[i + levelSize];
                 samples[i] = ai + t;
                 samples[i + levelSize] = ai - t;
             }
@@ -104,9 +101,9 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
             }
 
             Radix2Reorder(samples);
-            for (int levelSize = 1; levelSize < samples.Length; levelSize *= 2)
+            for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2)
             {
-                for (int k = 0; k < levelSize; k++)
+                for (var k = 0; k < levelSize; k++)
                 {
                     Radix2Step(samples, exponentSign, levelSize, k);
                 }
@@ -127,13 +124,18 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
             }
 
             Radix2Reorder(samples);
-            for (int levelSize = 1; levelSize < samples.Length; levelSize *= 2)
+            for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2)
             {
-                int size = levelSize;
-                Parallel.For(
-                    0,
-                    size,
-                    k => Radix2Step(samples, exponentSign, size, k));
+                var size = levelSize;
+                Parallel.ForEach(
+                    Partitioner.Create(0, size), 
+                    (range, loopState) =>
+                    {
+                        for (var k = range.Item1; k < range.Item2; k++)
+                        {
+                            Radix2Step(samples, exponentSign, size, k);
+                        }
+                    });
             }
         }
 
@@ -161,4 +163,4 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
             InverseScaleByOptions(options, samples);
         }
     }
-}
+}
\ No newline at end of file
diff --git a/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs b/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs
index 260031c6..2b079062 100644
--- a/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs
+++ b/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs
@@ -1,9 +1,7 @@
 // <copyright file="DiscreteHartleyTransform.Naive.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
-//
 // Copyright (c) 2009 Math.NET
-//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,10 +10,8 @@
 // copies of the Software, and to permit persons to whom the
 // Software is furnished to do so, subject to the following
 // conditions:
-//
 // The above copyright notice and this permission notice shall be
 // included in all copies or substantial portions of the Software.
-//
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
@@ -29,7 +25,8 @@
 namespace MathNet.Numerics.IntegralTransforms.Algorithms
 {
     using System;
-    using Threading;
+    using System.Collections.Concurrent;
+    using System.Threading.Tasks;
 
     /// <summary>
     /// Fast (FHT) Implementation of the Discrete Hartley Transform (DHT).
@@ -43,23 +40,25 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
         /// <returns>Corresponding frequency-space vector.</returns>
         internal static double[] Naive(double[] samples)
         {
-            double w0 = Constants.Pi2 / samples.Length;
+            var w0 = Constants.Pi2 / samples.Length;
             var spectrum = new double[samples.Length];
 
-            Parallel.For(
-                0,
-                samples.Length,
-                k =>
+            Parallel.ForEach(
+                Partitioner.Create(0, samples.Length), 
+                (range, loopState) =>
                 {
-                    double wk = w0 * k;
-                    double sum = 0.0;
-                    for (int n = 0; n < samples.Length; n++)
+                    for (var k = range.Item1; k < range.Item2; k++)
                     {
-                        double w = n * wk;
-                        sum += samples[n] * Constants.Sqrt2 * Math.Cos(w - Constants.PiOver4);
-                    }
+                        var wk = w0 * k;
+                        var sum = 0.0;
+                        for (var n = 0; n < samples.Length; n++)
+                        {
+                            var w = n * wk;
+                            sum += samples[n] * Constants.Sqrt2 * Math.Cos(w - Constants.PiOver4);
+                        }
 
-                    spectrum[k] = sum;
+                        spectrum[k] = sum;
+                    }
                 });
 
             return spectrum;
@@ -91,4 +90,4 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms
             return timeSpace;
         }
     }
-}
+}
\ No newline at end of file
diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
index 22b1c0df..10334373 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
@@ -1,9 +1,7 @@
 // <copyright file="DenseVector.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
-//
 // Copyright (c) 2009 Math.NET
-//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,10 +10,8 @@
 // copies of the Software, and to permit persons to whom the
 // Software is furnished to do so, subject to the following
 // conditions:
-//
 // The above copyright notice and this permission notice shall be
 // included in all copies or substantial portions of the Software.
-//
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
@@ -29,11 +25,12 @@
 namespace MathNet.Numerics.LinearAlgebra.Double
 {
     using System;
+    using System.Collections.Concurrent;
     using System.Collections.Generic;
     using System.Globalization;
+    using System.Threading.Tasks;
     using NumberTheory;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// A vector using dense storage.
@@ -52,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public DenseVector(int size)
             : base(size)
         {
-            Data = new double[size];
+            this.Data = new double[size];
         }
 
         /// <summary>
@@ -71,9 +68,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public DenseVector(int size, double value)
             : this(size)
         {
-            for (var index = 0; index < Data.Length; index++)
+            for (var index = 0; index < this.Data.Length; index++)
             {
-                Data[index] = value;
+                this.Data[index] = value;
             }
         }
 
@@ -90,11 +87,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             var vector = other as DenseVector;
             if (vector == null)
             {
-                Parallel.For(0, Count, index => this[index] = other[index]);
+                Parallel.ForEach(
+                    Partitioner.Create(0, this.Data.Length), 
+                    (range, loopState) =>
+                    {
+                        for (var index = range.Item1; index < range.Item2; index++)
+                        {
+                            this[index] = other[index];
+                        }
+                    });
             }
             else
             {
-                Buffer.BlockCopy(vector.Data, 0, Data, 0, Data.Length * Constants.SizeOfDouble);
+                Buffer.BlockCopy(vector.Data, 0, this.Data, 0, this.Data.Length * Constants.SizeOfDouble);
             }
         }
 
@@ -108,7 +113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public DenseVector(DenseVector other)
             : this(other.Count)
         {
-            Buffer.BlockCopy(other.Data, 0, Data, 0, Data.Length * Constants.SizeOfDouble);
+            Buffer.BlockCopy(other.Data, 0, this.Data, 0, this.Data.Length * Constants.SizeOfDouble);
         }
 
         /// <summary>
@@ -119,7 +124,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// changes to the vector will also change the array.</remarks>
         public DenseVector(double[] array) : base(array.Length)
         {
-            Data = array;
+            this.Data = array;
         }
 
         /// <summary>
@@ -174,10 +179,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <returns>This vector as a column matrix.</returns>
         public override Matrix ToColumnMatrix()
         {
-            var matrix = new DenseMatrix(Count, 1);
-            for (var i = 0; i < Data.Length; i++)
+            var matrix = new DenseMatrix(this.Count, 1);
+            for (var i = 0; i < this.Data.Length; i++)
             {
-                matrix[i, 0] = Data[i];
+                matrix[i, 0] = this.Data[i];
             }
 
             return matrix;
@@ -189,10 +194,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <returns>This vector as a row matrix.</returns>
         public override Matrix ToRowMatrix()
         {
-            var matrix = new DenseMatrix(1, Count);
-            for (var i = 0; i < Data.Length; i++)
+            var matrix = new DenseMatrix(1, this.Count);
+            for (var i = 0; i < this.Data.Length; i++)
             {
-                matrix[0, i] = Data[i];
+                matrix[0, i] = this.Data[i];
             }
 
             return matrix;
@@ -207,12 +212,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             get
             {
-                return Data[index];
+                return this.Data[index];
             }
 
             set
             {
-                Data[index] = value;
+                this.Data[index] = value;
             }
         }
 
@@ -268,11 +273,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("target");
             }
 
-            if (Count != target.Count)
+            if (this.Count != target.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "target");
             }
-            
+
             if (ReferenceEquals(this, target))
             {
                 return;
@@ -281,11 +286,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             var otherVector = target as DenseVector;
             if (otherVector == null)
             {
-                Parallel.For(0, Data.Length, index => target[index] = Data[index]);
+                Parallel.ForEach(
+                    Partitioner.Create(0, this.Data.Length), 
+                    (range, loopState) =>
+                    {
+                        for (var index = range.Item1; index < range.Item2; index++)
+                        {
+                            target[index] = this.Data[index];
+                        }
+                    });
             }
             else
             {
-                Buffer.BlockCopy(Data, 0, otherVector.Data, 0, Data.Length * Constants.SizeOfDouble);
+                Buffer.BlockCopy(this.Data, 0, otherVector.Data, 0, this.Data.Length * Constants.SizeOfDouble);
             }
         }
 
@@ -300,7 +313,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Parallel.For(0, Count, i => Data[i] += scalar);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Data.Length), 
+                (range, loopState) =>
+                {
+                    for (var index = range.Item1; index < range.Item2; index++)
+                    {
+                        this.Data[index] += scalar;
+                    }
+                });
         }
 
         /// <summary>
@@ -317,12 +338,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
 
-            CopyTo(result);
+            this.CopyTo(result);
             result.Add(scalar);
         }
 
@@ -339,7 +360,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
@@ -352,7 +373,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                Control.LinearAlgebraProvider.AddVectorToScaledVector(Data, 1.0, denseVector.Data);
+                Control.LinearAlgebraProvider.AddVectorToScaledVector(this.Data, 1.0, denseVector.Data);
             }
         }
 
@@ -372,12 +393,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
@@ -390,7 +411,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                CopyTo(result);
+                this.CopyTo(result);
                 result.Add(other);
             }
         }
@@ -453,7 +474,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Parallel.For(0, Count, i => Data[i] -= scalar);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Data.Length), 
+                (range, loopState) =>
+                {
+                    for (var index = range.Item1; index < range.Item2; index++)
+                    {
+                        this.Data[index] -= scalar;
+                    }
+                });
         }
 
         /// <summary>
@@ -470,12 +499,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
 
-            CopyTo(result);
+            this.CopyTo(result);
             result.Subtract(scalar);
         }
 
@@ -485,14 +514,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="other">The vector to subtract from this one.</param>
         /// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
         /// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
-         public override void Subtract(Vector other)
+        public override void Subtract(Vector other)
         {
             if (other == null)
             {
                 throw new ArgumentNullException("other");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
@@ -505,19 +534,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                Control.LinearAlgebraProvider.AddVectorToScaledVector(Data, -1.0, denseVector.Data);
+                Control.LinearAlgebraProvider.AddVectorToScaledVector(this.Data, -1.0, denseVector.Data);
             }
         }
 
-         /// <summary>
-         /// Subtracts another vector to this vector and stores the result into the result vector.
-         /// </summary>
-         /// <param name="other">The vector to subtract from this one.</param>
-         /// <param name="result">The vector to store the result of the subtraction.</param>
-         /// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
-         /// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
-         /// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
-         /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
+        /// <summary>
+        /// Subtracts another vector to this vector and stores the result into the result vector.
+        /// </summary>
+        /// <param name="other">The vector to subtract from this one.</param>
+        /// <param name="result">The vector to store the result of the subtraction.</param>
+        /// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
+        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
+        /// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
+        /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
         public override void Subtract(Vector other, Vector result)
         {
             if (result == null)
@@ -525,12 +554,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
@@ -543,7 +572,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                CopyTo(result);
+                this.CopyTo(result);
                 result.Subtract(other);
             }
         }
@@ -601,8 +630,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <remarks>Added as an alternative to the unary negation operator.</remarks>
         public override Vector Negate()
         {
-            var result = new DenseVector(Count);
-            Parallel.For(0, Count, i => result[i] = -Data[i]);
+            var result = new DenseVector(this.Count);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Data.Length), 
+                (range, loopState) =>
+                {
+                    for (var index = range.Item1; index < range.Item2; index++)
+                    {
+                        result[index] = -this.Data[index];
+                    }
+                });
+
             return result;
         }
 
@@ -617,7 +655,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Control.LinearAlgebraProvider.ScaleArray(scalar, Data);
+            Control.LinearAlgebraProvider.ScaleArray(scalar, this.Data);
         }
 
         /// <summary>
@@ -634,22 +672,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-
             var denseVector = other as DenseVector;
 
             if (denseVector == null)
             {
                 return base.DotProduct(other);
             }
-            else
-            {
-                return Control.LinearAlgebraProvider.DotProduct(this.Data, denseVector.Data);
-            }
+            
+            return Control.LinearAlgebraProvider.DotProduct(this.Data, denseVector.Data);
         }
 
         /// <summary>
@@ -746,9 +781,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public override double Norm()
         {
             var sum = 0.0;
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < this.Data.Length; i++)
             {
-                sum = SpecialFunctions.Hypotenuse(sum, Data[i]);
+                sum = SpecialFunctions.Hypotenuse(sum, this.Data[i]);
             }
 
             return sum;
@@ -762,16 +797,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             var sum = 0.0;
             var syncLock = new object();
-            
-            Parallel.For(
-                0, 
-                Count,
-                () => 0.0,
-                (index, localData) =>
+
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                () => 0.0, 
+                (range, loop, localData) =>
                 {
-                    localData += Math.Abs(Data[index]);
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        localData += Math.Abs(this.Data[i]);
+                    }
+
                     return localData;
-                },
+                }, 
                 localResult =>
                 {
                     lock (syncLock)
@@ -797,26 +835,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             if (1 == p)
             {
-                return Norm1();
+                return this.Norm1();
             }
 
             if (2 == p)
             {
-                return Norm();
+                return this.Norm();
             }
 
             var sum = 0.0;
             var syncLock = new object();
-           
-            Parallel.For(
-                0, 
-                Count,
-                () => 0.0,
-                (index, localData) =>
+
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                () => 0.0, 
+                (range, loop, localData) =>
                 {
-                    localData += Math.Pow(Math.Abs(Data[index]), p);
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        localData += Math.Pow(Math.Abs(this.Data[i]), p);
+                    }
+
                     return localData;
-                },
+                }, 
                 localResult =>
                 {
                     lock (syncLock)
@@ -838,15 +879,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             var syncLock = new object();
 
-            Parallel.For(
-                0, 
-                Count,
-                () => 0.0,
-                (index, localData) =>
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                () => 0.0, 
+                (range, loop, localData) =>
                 {
-                    localData = Math.Max(localData, Math.Abs(Data[index]));
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        localData = Math.Max(localData, Math.Abs(this.Data[i]));
+                    }
+
                     return localData;
-                },
+                }, 
                 localResult =>
                 {
                     lock (syncLock)
@@ -940,7 +984,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             // parsing
             var data = new double[(tokens.Count + 1) >> 1];
-            for (int i = 0; i < data.Length; i++)
+            for (var i = 0; i < data.Length; i++)
             {
                 if (token == null || token.Value == textInfo.ListSeparator)
                 {
diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs
index e219cd22..49ad3c71 100644
--- a/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs
+++ b/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs
@@ -1,9 +1,7 @@
 // <copyright file="Matrix.Arithmetic.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
-//
 // Copyright (c) 2009 Math.NET
-//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,10 +10,8 @@
 // copies of the Software, and to permit persons to whom the
 // Software is furnished to do so, subject to the following
 // conditions:
-//
 // The above copyright notice and this permission notice shall be
 // included in all copies or substantial portions of the Software.
-//
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
@@ -29,9 +25,9 @@
 namespace MathNet.Numerics.LinearAlgebra.Double
 {
     using System;
-    using System.Text;
+    using System.Collections.Concurrent;
+    using System.Threading.Tasks;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// Defines the base class for <c>Matrix</c> classes.
@@ -51,19 +47,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
+            if (other.RowCount != this.RowCount || other.ColumnCount != this.ColumnCount)
             {
                 throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions);
             }
 
-            Parallel.For(
-                0, 
-                RowCount, 
-                i =>
+            Parallel.ForEach(
+                Partitioner.Create(0, this.RowCount), 
+                (range, loopState) =>
                 {
-                    for (int j = 0; j < ColumnCount; j++)
+                    for (var i = range.Item1; i < range.Item2; i++)
                     {
-                        At(i, j, At(i, j) + other.At(i, j));
+                        for (var j = 0; j < this.ColumnCount; j++)
+                        {
+                            this.At(i, j, this.At(i, j) + other.At(i, j));
+                        }
                     }
                 });
         }
@@ -81,19 +79,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
+            if (other.RowCount != this.RowCount || other.ColumnCount != this.ColumnCount)
             {
                 throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions);
             }
 
-            Parallel.For(
-                0,
-                RowCount,
-                i =>
+            Parallel.ForEach(
+                Partitioner.Create(0, this.RowCount), 
+                (range, loopState) =>
                 {
-                    for (int j = 0; j < ColumnCount; j++)
+                    for (var i = range.Item1; i < range.Item2; i++)
                     {
-                        At(i, j, At(i, j) - other.At(i, j));
+                        for (var j = 0; j < this.ColumnCount; j++)
+                        {
+                            this.At(i, j, this.At(i, j) - other.At(i, j));
+                        }
                     }
                 });
         }
@@ -104,19 +104,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="scalar">The scalar to multiply with.</param>
         public virtual void Multiply(double scalar)
         {
-            if (Precision.AlmostEqualInDecimalPlaces(1.0, scalar, 15))
+            if (1.0.AlmostEqualInDecimalPlaces(scalar, 15))
             {
                 return;
             }
 
-            Parallel.For(
-                0,
-                RowCount,
-                i =>
+            Parallel.ForEach(
+                Partitioner.Create(0, this.RowCount), 
+                (range, loopState) =>
                 {
-                    for (int j = 0; j < ColumnCount; j++)
+                    for (var i = range.Item1; i < range.Item2; i++)
                     {
-                        At(i, j, At(i, j) * scalar);
+                        for (var j = 0; j < this.ColumnCount; j++)
+                        {
+                            this.At(i, j, this.At(i, j) * scalar);
+                        }
                     }
                 });
         }
@@ -135,17 +137,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (result.RowCount != RowCount)
+            if (result.RowCount != this.RowCount)
             {
-                throw new ArgumentException("result", Resources.ArgumentMatrixSameRowDimension);
+                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension, "result");
             }
 
-            if (result.ColumnCount != ColumnCount)
+            if (result.ColumnCount != this.ColumnCount)
             {
-                throw new ArgumentException("result", Resources.ArgumentMatrixSameColumnDimension);
+                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "result");
             }
 
-            CopyTo(result);
+            this.CopyTo(result);
             result.Multiply(scalar);
         }
 
@@ -158,7 +160,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentException">If <c>this.ColumnCount != rightSide.Count</c>.</exception>
         public virtual Vector Multiply(Vector rightSide)
         {
-            Vector ret = CreateVector(RowCount);
+            var ret = this.CreateVector(this.RowCount);
             Multiply(rightSide, ret);
             return ret;
         }
@@ -179,9 +181,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("rightSide");
             }
 
-            if (ColumnCount != rightSide.Count)
+            if (this.ColumnCount != rightSide.Count)
             {
-                throw new ArgumentException("rightSide", Resources.ArgumentMatrixDimensions);
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions, "rightSide");
             }
 
             if (result == null)
@@ -189,31 +191,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (RowCount != result.Count)
+            if (this.RowCount != result.Count)
             {
-                throw new ArgumentException("result", Resources.ArgumentMatrixDimensions);
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
             }
 
             if (ReferenceEquals(rightSide, result))
             {
-                Vector tmp = result.CreateVector(result.Count);
+                var tmp = result.CreateVector(result.Count);
                 Multiply(rightSide, tmp);
                 tmp.CopyTo(result);
             }
             else
             {
-                Parallel.For(
-                    0,
-                    RowCount,
-                    i =>
+                Parallel.ForEach(
+                    Partitioner.Create(0, this.RowCount), 
+                    (range, loopState) =>
                     {
-                        double s = 0;
-                        for (int j = 0; j != ColumnCount; j++)
+                        for (var i = range.Item1; i < range.Item2; i++)
                         {
-                            s += At(i, j) * rightSide[j];
-                        }
+                            double s = 0;
+                            for (var j = 0; j != this.ColumnCount; j++)
+                            {
+                                s += this.At(i, j) * rightSide[j];
+                            }
 
-                        result[i] = s;
+                            result[i] = s;
+                        }
                     });
             }
         }
@@ -227,8 +231,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentException">If <strong>this.RowCount != leftSide.Count</strong>.</exception>
         public virtual Vector LeftMultiply(Vector leftSide)
         {
-            Vector ret = CreateVector(ColumnCount);
-            LeftMultiply(leftSide, ret);
+            var ret = this.CreateVector(this.ColumnCount);
+            this.LeftMultiply(leftSide, ret);
             return ret;
         }
 
@@ -248,9 +252,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("leftSide");
             }
 
-            if (RowCount != leftSide.Count)
+            if (this.RowCount != leftSide.Count)
             {
-                throw new ArgumentException("leftSide", Resources.ArgumentMatrixDimensions);
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions, "leftSide");
             }
 
             if (result == null)
@@ -258,31 +262,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (ColumnCount != result.Count)
+            if (this.ColumnCount != result.Count)
             {
-                throw new ArgumentException("result", Resources.ArgumentMatrixDimensions);
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
             }
 
             if (ReferenceEquals(leftSide, result))
             {
-                Vector tmp = result.CreateVector(result.Count);
-                LeftMultiply(leftSide, tmp);
+                var tmp = result.CreateVector(result.Count);
+                this.LeftMultiply(leftSide, tmp);
                 tmp.CopyTo(result);
             }
             else
             {
-                Parallel.For(
-                    0,
-                    ColumnCount,
-                    j =>
+                Parallel.ForEach(
+                    Partitioner.Create(0, this.ColumnCount), 
+                    (range, loopState) =>
                     {
-                        double s = 0;
-                        for (int i = 0; i != leftSide.Count; i++)
+                        for (var j = range.Item1; j < range.Item2; j++)
                         {
-                            s += leftSide[i] * At(i, j);
-                        }
+                            double s = 0;
+                            for (var i = 0; i != leftSide.Count; i++)
+                            {
+                                s += leftSide[i] * this.At(i, j);
+                            }
 
-                        result[j] = s;
+                            result[j] = s;
+                        }
                     });
             }
         }
@@ -308,38 +314,40 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (ColumnCount != other.RowCount)
+            if (this.ColumnCount != other.RowCount)
             {
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
-            if (result.RowCount != RowCount || result.ColumnCount != other.ColumnCount)
+            if (result.RowCount != this.RowCount || result.ColumnCount != other.ColumnCount)
             {
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
             if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
             {
-                Matrix tmp = result.CreateMatrix(result.RowCount, result.ColumnCount);
+                var tmp = result.CreateMatrix(result.RowCount, result.ColumnCount);
                 Multiply(other, tmp);
                 tmp.CopyTo(result);
             }
             else
             {
-                Parallel.For(
-                    0,
-                    this.RowCount,
-                    j =>
+                Parallel.ForEach(
+                    Partitioner.Create(0, this.RowCount), 
+                    (range, loopState) =>
                     {
-                        for (int i = 0; i != other.ColumnCount; i++)
+                        for (var j = range.Item1; j < range.Item2; j++)
                         {
-                            double s = 0;
-                            for (int l = 0; l < this.ColumnCount; l++)
+                            for (var i = 0; i != other.ColumnCount; i++)
                             {
-                                s += this.At(j, l) * other.At(l, i);
-                            }
+                                double s = 0;
+                                for (var l = 0; l < this.ColumnCount; l++)
+                                {
+                                    s += this.At(j, l) * other.At(l, i);
+                                }
 
-                            result.At(j, i, s);
+                                result.At(j, i, s);
+                            }
                         }
                     });
             }
@@ -359,12 +367,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (ColumnCount != other.RowCount)
+            if (this.ColumnCount != other.RowCount)
             {
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
-            Matrix result = CreateMatrix(RowCount, other.ColumnCount);
+            var result = this.CreateMatrix(this.RowCount, other.ColumnCount);
             Multiply(other, result);
             return result;
         }
@@ -374,7 +382,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>        
         public virtual void Negate()
         {
-            Multiply(-1);
+            this.Multiply(-1);
         }
 
         /// <summary>
@@ -390,12 +398,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
+            if (result.RowCount != this.RowCount || result.ColumnCount != this.ColumnCount)
             {
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
-            CopyTo(result);
+            this.CopyTo(result);
             result.Negate();
         }
 
@@ -427,7 +435,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions);
             }
 
-            Matrix ret = leftSide.Clone();
+            var ret = leftSide.Clone();
             ret.Add(rightSide);
             return ret;
         }
@@ -476,7 +484,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions);
             }
 
-            Matrix ret = leftSide.Clone();
+            var ret = leftSide.Clone();
             ret.Subtract(rightSide);
             return ret;
         }
@@ -494,7 +502,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("rightSide");
             }
 
-            Matrix ret = rightSide.Clone();
+            var ret = rightSide.Clone();
             ret.Negate();
             return ret;
         }
@@ -513,7 +521,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("leftSide");
             }
 
-            Matrix ret = leftSide.Clone();
+            var ret = leftSide.Clone();
             ret.Multiply(rightSide);
             return ret;
         }
@@ -532,7 +540,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("rightSide");
             }
 
-            Matrix ret = rightSide.Clone();
+            var ret = rightSide.Clone();
             ret.Multiply(leftSide);
             return ret;
         }
diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs
index b53641f5..e8cf124c 100644
--- a/src/Numerics/LinearAlgebra/Double/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs
@@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
     using System;
     using System.Text;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// Defines the base class for <c>Matrix</c> classes.
@@ -270,7 +269,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
         /// or greater than or equal to the number of rows.</exception>
-        /// <exception cref="NotConformableException">If <b>this.Columns != result.Count</b>.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If <b>this.Columns != result.Count</b>.</exception>
         public virtual void GetRow(int index, Vector result)
         {
             GetRow(index, 0, ColumnCount, result);
@@ -312,7 +311,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> + <paramref name="length"/>  
         /// is greater than or equal to the number of rows.</exception>
         /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
-        /// <exception cref="NotConformableException">If <strong>result.Count &lt; length</strong>.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If <strong>result.Count &lt; length</strong>.</exception>
         public virtual void GetRow(int rowIndex, int columnIndex, int length, Vector result)
         {
             if (result == null)
@@ -324,6 +323,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             {
                 throw new ArgumentOutOfRangeException("rowIndex");
             }
+
             if (columnIndex >= ColumnCount || columnIndex < 0)
             {
                 throw new ArgumentOutOfRangeException("columnIndex");
@@ -341,7 +341,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             if (result.Count < length)
             {
-                throw new ArgumentException("result", Resources.ArgumentVectorsSameLength);
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
 
             for (int i = columnIndex, j = 0; i < columnIndex + length; i++, j++)
@@ -372,7 +372,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentNullException">If the result <see cref="Vector"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
         /// or greater than or equal to the number of columns.</exception>
-        /// <exception cref="NotConformableException">If <b>this.Rows != result.Count</b>.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If <b>this.Rows != result.Count</b>.</exception>
         public virtual void GetColumn(int index, Vector result)
         {
             GetColumn(index, 0, RowCount, result);
@@ -415,7 +415,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> + <paramref name="length"/>  
         /// is greater than or equal to the number of rows.</exception>
         /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
-        /// <exception cref="NotConformableException">If <strong>result.Count &lt; length</strong>.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If <strong>result.Count &lt; length</strong>.</exception>
         public virtual void GetColumn(int columnIndex, int rowIndex, int length, Vector result)
         {
             if (result == null)
@@ -427,6 +427,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             {
                 throw new ArgumentOutOfRangeException("columnIndex");
             }
+
             if (rowIndex >= RowCount || rowIndex < 0)
             {
                 throw new ArgumentOutOfRangeException("rowIndex");
@@ -444,7 +445,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             if (result.Count < length)
             {
-                throw new ArgumentException("result", Resources.ArgumentVectorsSameLength);
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
 
             for (int i = rowIndex, j = 0; i < rowIndex + length; i++, j++)
diff --git a/src/Numerics/LinearAlgebra/Double/Vector.cs b/src/Numerics/LinearAlgebra/Double/Vector.cs
index 88dabe95..2dc077b8 100644
--- a/src/Numerics/LinearAlgebra/Double/Vector.cs
+++ b/src/Numerics/LinearAlgebra/Double/Vector.cs
@@ -1,9 +1,7 @@
 // <copyright file="Vector.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
+// Copyright (c) 2009-2010 Math.NET
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,10 +10,8 @@
 // copies of the Software, and to permit persons to whom the
 // Software is furnished to do so, subject to the following
 // conditions:
-//
 // The above copyright notice and this permission notice shall be
 // included in all copies or substantial portions of the Software.
-//
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
@@ -30,11 +26,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 {
     using System;
     using System.Collections;
+    using System.Collections.Concurrent;
     using System.Collections.Generic;
     using System.Text;
-
+    using System.Threading.Tasks;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// Defines the base class for <c>Vector</c> classes.
@@ -46,8 +42,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 #if SILVERLIGHT
     IFormattable, IEnumerable<double>, IEquatable<Vector>
 #else
-    IFormattable, IEnumerable<double>, IEquatable<Vector>, ICloneable
-#endif    
+        IFormattable, IEnumerable<double>, IEquatable<Vector>, ICloneable
+#endif
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="Vector"/> class. 
@@ -66,7 +62,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentException(Resources.ArgumentMustBePositive, "size");
             }
 
-            Count = size;
+            this.Count = size;
         }
 
         /// <summary>
@@ -117,10 +113,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public abstract Vector CreateVector(int size);
 
         #region Elementary operations
+
         /// <summary>
         /// Adds a scalar to each element of the vector.
         /// </summary>
-        /// <param name="scalar">The scalar to add.</param>
+        /// <param name="scalar">
+        /// The scalar to add.
+        /// </param>
         public virtual void Add(double scalar)
         {
             if (scalar == 0.0)
@@ -128,16 +127,32 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Parallel.For(0, Count, i => this[i] += scalar);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        this[i] += scalar;
+                    }
+                });
         }
 
         /// <summary>
         /// Adds a scalar to each element of the vector and stores the result in the result vector.
         /// </summary>
-        /// <param name="scalar">The scalar to add.</param>
-        /// <param name="result">The vector to store the result of the addition.</param>
-        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
+        /// <param name="scalar">
+        /// The scalar to add.
+        /// </param>
+        /// <param name="result">
+        /// The vector to store the result of the addition.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the result vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="result"/> are not the same size.
+        /// </exception>
         public virtual void Add(double scalar, Vector result)
         {
             if (result == null)
@@ -145,14 +160,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
 
             if (!ReferenceEquals(this, result))
             {
-                CopyTo(result);
+                this.CopyTo(result);
             }
 
             result.Add(scalar);
@@ -161,8 +176,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Returns this vector.
         /// </summary>
-        /// <returns>This vector.</returns>
-        /// <remarks>Added as an alternative to the unary addition operator.</remarks>
+        /// <returns>
+        /// This vector.
+        /// </returns>
+        /// <remarks>
+        /// Added as an alternative to the unary addition operator.
+        /// </remarks>
         public virtual Vector Plus()
         {
             return this;
@@ -171,9 +190,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Adds another vector to this vector.
         /// </summary>
-        /// <param name="other">The vector to add to this one.</param>
-        /// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception> 
-        /// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
+        /// <param name="other">
+        /// The vector to add to this one.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the other vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="other"/> are not the same size.
+        /// </exception>
         public virtual void Add(Vector other)
         {
             if (other == null)
@@ -181,23 +206,43 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            Parallel.For(0, Count, i => this[i] += other[i]);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        this[i] += other[i];
+                    }
+                });
         }
 
         /// <summary>
         /// Adds another vector to this vector and stores the result into the result vector.
         /// </summary>
-        /// <param name="other">The vector to add to this one.</param>
-        /// <param name="result">The vector to store the result of the addition.</param>
-        /// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception> 
-        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception> 
-        /// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
-        /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
+        /// <param name="other">
+        /// The vector to add to this one.
+        /// </param>
+        /// <param name="result">
+        /// The vector to store the result of the addition.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the other vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentNullException">
+        /// If the result vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="other"/> are not the same size.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="result"/> are not the same size.
+        /// </exception>
         public virtual void Add(Vector other, Vector result)
         {
             if (result == null)
@@ -205,7 +250,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
@@ -218,7 +263,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                CopyTo(result);
+                this.CopyTo(result);
                 result.Add(other);
             }
         }
@@ -226,7 +271,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Subtracts a scalar from each element of the vector.
         /// </summary>
-        /// <param name="scalar">The scalar to subtract.</param>
+        /// <param name="scalar">
+        /// The scalar to subtract.
+        /// </param>
         public virtual void Subtract(double scalar)
         {
             if (scalar == 0.0)
@@ -234,16 +281,32 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Parallel.For(0, Count, i => this[i] -= scalar);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        this[i] -= scalar;
+                    }
+                });
         }
 
         /// <summary>
-        ///  Subtracts a scalar from each element of the vector and stores the result in the result vector.
+        /// Subtracts a scalar from each element of the vector and stores the result in the result vector.
         /// </summary>
-        /// <param name="scalar">The scalar to subtract.</param>
-        /// <param name="result">The vector to store the result of the subtraction.</param>
-        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
+        /// <param name="scalar">
+        /// The scalar to subtract.
+        /// </param>
+        /// <param name="result">
+        /// The vector to store the result of the subtraction.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the result vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="result"/> are not the same size.
+        /// </exception>
         public virtual void Subtract(double scalar, Vector result)
         {
             if (result == null)
@@ -251,14 +314,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
 
             if (!ReferenceEquals(this, result))
             {
-                CopyTo(result);
+                this.CopyTo(result);
             }
 
             result.Subtract(scalar);
@@ -267,8 +330,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Returns a negated vector.
         /// </summary>
-        /// <returns>The negated vector.</returns>
-        /// <remarks>Added as an alternative to the unary negation operator.</remarks>
+        /// <returns>
+        /// The negated vector.
+        /// </returns>
+        /// <remarks>
+        /// Added as an alternative to the unary negation operator.
+        /// </remarks>
         public virtual Vector Negate()
         {
             return this * -1;
@@ -277,9 +344,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Subtracts another vector from this vector.
         /// </summary>
-        /// <param name="other">The vector to subtract from this one.</param>
-        /// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception> 
-        /// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
+        /// <param name="other">
+        /// The vector to subtract from this one.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the other vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="other"/> are not the same size.
+        /// </exception>
         public virtual void Subtract(Vector other)
         {
             if (other == null)
@@ -287,23 +360,43 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            Parallel.For(0, Count, i => this[i] -= other[i]);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        this[i] -= other[i];
+                    }
+                });
         }
 
         /// <summary>
         /// Subtracts another vector to this vector and stores the result into the result vector.
         /// </summary>
-        /// <param name="other">The vector to subtract from this one.</param>
-        /// <param name="result">The vector to store the result of the subtraction.</param>
-        /// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception> 
-        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception> 
-        /// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
-        /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
+        /// <param name="other">
+        /// The vector to subtract from this one.
+        /// </param>
+        /// <param name="result">
+        /// The vector to store the result of the subtraction.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the other vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentNullException">
+        /// If the result vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="other"/> are not the same size.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="result"/> are not the same size.
+        /// </exception>
         public virtual void Subtract(Vector other, Vector result)
         {
             if (result == null)
@@ -311,7 +404,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
@@ -324,7 +417,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                CopyTo(result);
+                this.CopyTo(result);
                 result.Subtract(other);
             }
         }
@@ -332,7 +425,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Multiplies a scalar to each element of the vector.
         /// </summary>
-        /// <param name="scalar">The scalar to multiply.</param>
+        /// <param name="scalar">
+        /// The scalar to multiply.
+        /// </param>
         public virtual void Multiply(double scalar)
         {
             if (scalar == 1.0)
@@ -340,16 +435,32 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Parallel.For(0, Count, index => this[index] *= scalar);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                (range, loopState) =>
+                {
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        this[i] *= scalar;
+                    }
+                });
         }
 
         /// <summary>
         /// Multiplies a scalar to each element of the vector and stores the result in the result vector.
         /// </summary>
-        /// <param name="scalar">The scalar to multiply.</param>
-        /// <param name="result">The vector to store the result of the multiplication.</param>
-        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
+        /// <param name="scalar">
+        /// The scalar to multiply.
+        /// </param>
+        /// <param name="result">
+        /// The vector to store the result of the multiplication.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the result vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="result"/> are not the same size.
+        /// </exception>
         public virtual void Multiply(double scalar, Vector result)
         {
             if (result == null)
@@ -357,14 +468,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
-        
+
             if (!ReferenceEquals(this, result))
             {
-                CopyTo(result);
+                this.CopyTo(result);
             }
 
             result.Multiply(scalar);
@@ -373,10 +484,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Computes the dot product between this vector and another vector.
         /// </summary>
-        /// <param name="other">The other vector to add.</param>
-        /// <returns>The result of the addition.</returns>
-        /// <exception cref="ArgumentException">If <paramref name="other"/> is not of the same size.</exception>
-        /// <exception cref="ArgumentNullException">If <paramref name="other"/> is <see langword="null" />.</exception>
+        /// <param name="other">
+        /// The other vector to add.
+        /// </param>
+        /// <returns>
+        /// The result of the addition.
+        /// </returns>
+        /// <exception cref="ArgumentException">
+        /// If <paramref name="other"/> is not of the same size.
+        /// </exception>
+        /// <exception cref="ArgumentNullException">
+        /// If <paramref name="other"/> is <see langword="null"/>.
+        /// </exception>
         public virtual double DotProduct(Vector other)
         {
             if (other == null)
@@ -384,13 +503,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("other");
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            double dot = 0.0;
-            for (int i = 0; i < Count; i++)
+            var dot = 0.0;
+            for (var i = 0; i < this.Count; i++)
             {
                 dot += this[i] * other[i];
             }
@@ -401,7 +520,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Divides each element of the vector by a scalar.
         /// </summary>
-        /// <param name="scalar">The scalar to divide with.</param>
+        /// <param name="scalar">
+        /// The scalar to divide with.
+        /// </param>
         public virtual void Divide(double scalar)
         {
             if (scalar == 1.0)
@@ -409,16 +530,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Multiply(1.0 / scalar);
+            this.Multiply(1.0 / scalar);
         }
 
         /// <summary>
-        ///  Divides each element of the vector by a scalar and stores the result in the result vector.
+        /// Divides each element of the vector by a scalar and stores the result in the result vector.
         /// </summary>
-        /// <param name="scalar">The scalar to divide with.</param>
-        /// <param name="result">The vector to store the result of the division.</param>
-        /// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
+        /// <param name="scalar">
+        /// The scalar to divide with.
+        /// </param>
+        /// <param name="result">
+        /// The vector to store the result of the division.
+        /// </param>
+        /// <exception cref="ArgumentNullException">
+        /// If the result vector is <see langword="null"/>.
+        /// </exception>
+        /// <exception cref="ArgumentException">
+        /// If this vector and <paramref name="result"/> are not the same size.
+        /// </exception>
         public virtual void Divide(double scalar, Vector result)
         {
             if (result == null)
@@ -426,21 +555,23 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("result");
             }
 
-            if (Count != result.Count)
+            if (this.Count != result.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
             }
 
             if (!ReferenceEquals(this, result))
             {
-                CopyTo(result);
+                this.CopyTo(result);
             }
 
             result.Multiply(1.0 / scalar);
         }
+
         #endregion
 
         #region Arithmetic Operator Overloading
+
         /// <summary>
         /// Returns a <strong>Vector</strong> containing the same values of rightSide. 
         /// </summary>
@@ -631,7 +762,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public virtual double Norm()
         {
-            return NormP(2);
+            return this.NormP(2);
         }
 
         /// <summary>
@@ -642,7 +773,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public virtual double SquaredNorm()
         {
-            var norm = Norm();
+            var norm = this.Norm();
             return norm * norm;
         }
 
@@ -654,14 +785,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public virtual double Norm1()
         {
-            return NormP(1);
+            return this.NormP(1);
         }
 
         /// <summary>
         /// Computes the p-Norm.
         /// </summary>
-        /// <param name="p">The p value.</param>
-        /// <returns>Scalar ret = (sum(abs(this[i])^p))^(1/p)</returns>
+        /// <param name="p">
+        /// The p value.
+        /// </param>
+        /// <returns>
+        /// Scalar ret = (sum(abs(this[i])^p))^(1/p)
+        /// </returns>
         public virtual double NormP(int p)
         {
             if (1 > p)
@@ -671,15 +806,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             var sum = 0.0;
             var syncLock = new object();
-            Parallel.For(
-                0, 
-                Count, 
-                () => 0.0,
-                (index, localData) =>
+
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                () => 0.0, 
+                (range, loop, localData) =>
                 {
-                    localData += Math.Pow(Math.Abs(this[index]), p);
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        localData += Math.Pow(Math.Abs(this[i]), p);
+                    }
+
                     return localData;
-                },
+                }, 
                 localResult =>
                 {
                     lock (syncLock)
@@ -701,15 +840,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             var max = 0.0;
             var syncLock = new object();
-            Parallel.For(
-                0, 
-                Count,
-                () => 0.0,
-                (index, localData) =>
+
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                () => 0.0, 
+                (range, loop, localData) =>
                 {
-                    localData = Math.Max(localData, Math.Abs(this[index]));
+                    for (var i = range.Item1; i < range.Item2; i++)
+                    {
+                        localData = Math.Max(localData, Math.Abs(this[i]));
+                    }
+
                     return localData;
-                },
+                }, 
                 localResult =>
                 {
                     lock (syncLock)
@@ -724,11 +867,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Normalizes this vector to a unit vector with respect to the Eucliden 2-Norm.
         /// </summary>
-        /// <returns>This vector normalized to a unit vector with respect to the Eucliden 2-Norm.</returns>
+        /// <returns>
+        /// This vector normalized to a unit vector with respect to the Eucliden 2-Norm.
+        /// </returns>
         public virtual Vector Normalize()
         {
-            var norm = Norm();
-            var clone = Clone();
+            var norm = this.Norm();
+            var clone = this.Clone();
             if (norm == 0.0)
             {
                 return clone;
@@ -741,6 +886,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         #endregion
 
         #region Copying and Conversion
+
         /// <summary>
         /// Returns a deep-copy clone of the vector.
         /// </summary>
@@ -749,8 +895,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public Vector Clone()
         {
-            var retrunVector = CreateVector(Count);
-            CopyTo(retrunVector);
+            var retrunVector = this.CreateVector(this.Count);
+            this.CopyTo(retrunVector);
             return retrunVector;
         }
 
@@ -773,7 +919,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("target");
             }
 
-            if (Count != target.Count)
+            if (this.Count != target.Count)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "target");
             }
@@ -783,7 +929,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            Parallel.For(0, Count, index => target[index] = this[index]);
+            Parallel.ForEach(
+                Partitioner.Create(0, this.Count), 
+                (range, loopState) =>
+                {
+                    for (var index = range.Item1; index < range.Item2; index++)
+                    {
+                        target[index] = this[index];
+                    }
+                });
         }
 
         /// <summary>
@@ -808,12 +962,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("destination");
             }
 
-            if (offset >= Count)
+            if (offset >= this.Count)
             {
                 throw new ArgumentOutOfRangeException("offset");
             }
 
-            if (offset + count > Count)
+            if (offset + count > this.Count)
             {
                 throw new ArgumentOutOfRangeException("count");
             }
@@ -831,22 +985,41 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             if (ReferenceEquals(this, destination))
             {
                 var tmpVector = destination.CreateVector(destination.Count);
-                CopyTo(tmpVector);
-                Parallel.For(0, count, index => destination[destinationOffset + index] = this[offset + index]);
+                this.CopyTo(tmpVector);
+
+                Parallel.ForEach(
+                    Partitioner.Create(0, count), 
+                    (range, loopState) =>
+                    {
+                        for (var index = range.Item1; index < range.Item2; index++)
+                        {
+                            destination[destinationOffset + index] = tmpVector[offset + index];
+                        }
+                    });
             }
             else
             {
-               Parallel.For(0, count, index => destination[destinationOffset + index] = this[offset + index]);
+                Parallel.ForEach(
+                    Partitioner.Create(0, count), 
+                    (range, loopState) =>
+                    {
+                        for (var index = range.Item1; index < range.Item2; index++)
+                        {
+                            destination[destinationOffset + index] = this[offset + index];
+                        }
+                    });
             }
         }
 
         /// <summary>
         /// Returns the data contained in the vector as an array.
         /// </summary>
-        /// <returns>The vector's data as an array.</returns>        
+        /// <returns>
+        /// The vector's data as an array.
+        /// </returns>
         public virtual double[] ToArray()
         {
-            var ret = new double[Count];
+            var ret = new double[this.Count];
             for (var i = 0; i < ret.Length; i++)
             {
                 ret[i] = this[i];
@@ -858,11 +1031,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Create a matrix based on this vector in column form (one single column).
         /// </summary>
-        /// <returns>This vector as a column matrix.</returns>
+        /// <returns>
+        /// This vector as a column matrix.
+        /// </returns>
         public virtual Matrix ToColumnMatrix()
         {
-            var matrix = CreateMatrix(Count, 1);
-            for (var i = 0; i < Count; i++)
+            var matrix = this.CreateMatrix(this.Count, 1);
+            for (var i = 0; i < this.Count; i++)
             {
                 matrix[i, 0] = this[i];
             }
@@ -873,11 +1048,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Create a matrix based on this vector in row form (one single row).
         /// </summary>
-        /// <returns>This vector as a row matrix.</returns>
+        /// <returns>
+        /// This vector as a row matrix.
+        /// </returns>
         public virtual Matrix ToRowMatrix()
         {
-            var matrix = CreateMatrix(1, Count);
-            for (var i = 0; i < Count; i++)
+            var matrix = this.CreateMatrix(1, this.Count);
+            for (var i = 0; i < this.Count; i++)
             {
                 matrix[0, i] = this[i];
             }
@@ -890,6 +1067,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         #region Implemented Interfaces
 
 #if !SILVERLIGHT
+
         #region ICloneable
 
         /// <summary>
@@ -900,10 +1078,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         object ICloneable.Clone()
         {
-            return Clone();
+            return this.Clone();
         }
 
         #endregion
+
 #endif
 
         #region IEnumerable
@@ -916,7 +1095,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         IEnumerator IEnumerable.GetEnumerator()
         {
-            return GetEnumerator();
+            return this.GetEnumerator();
         }
 
         #endregion
@@ -931,29 +1110,36 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public virtual IEnumerator<double> GetEnumerator()
         {
-            for (var index = 0; index < Count; index++)
+            for (var index = 0; index < this.Count; index++)
             {
                 yield return this[index];
             }
         }
 
         #endregion
+
         /// <summary>
         /// Returns an <see cref="IEnumerator{T}"/> that contains the position and value of the element.
         /// </summary>
-        /// <returns>An <see cref="IEnumerator{T}"/> over this vector that contains the position and value of each
-        /// non-zero element.</returns>
-        /// <remarks>The enumerator returns a <seealso cref="KeyValuePair{T,K}"/> with the key being the element index and the value 
+        /// <returns>
+        /// An <see cref="IEnumerator{T}"/> over this vector that contains the position and value of each
+        /// non-zero element.
+        /// </returns>
+        /// <remarks>
+        /// The enumerator returns a 
+        /// <seealso cref="KeyValuePair{T,K}"/>
+        /// with the key being the element index and the value 
         /// being the value of the element at that index. For sparse vectors, the enumerator will exclude all elements
-        /// with a zero value.</remarks>
+        /// with a zero value.
+        /// </remarks>
         public virtual IEnumerable<KeyValuePair<int, double>> GetIndexedEnumerator()
         {
-            for (int i = 0; i < Count; i++)
+            for (var i = 0; i < this.Count; i++)
             {
                 yield return new KeyValuePair<int, double>(i, this[i]);
             }
         }
- 
+
         #region IEquatable<Vector>
 
         /// <summary>
@@ -973,7 +1159,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return false;
             }
 
-            if (Count != other.Count)
+            if (this.Count != other.Count)
             {
                 return false;
             }
@@ -985,7 +1171,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             // If all else fails, perform element wise comparison.
-            for (var index = 0; index < Count; index++)
+            for (var index = 0; index < this.Count; index++)
             {
                 if (!this[index].AlmostEqual(other[index]))
                 {
@@ -1011,7 +1197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public string ToString(IFormatProvider formatProvider)
         {
-            return ToString(null, formatProvider);
+            return this.ToString(null, formatProvider);
         }
 
         /// <summary>
@@ -1029,10 +1215,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public string ToString(string format, IFormatProvider formatProvider)
         {
             var stringBuilder = new StringBuilder();
-            for (var index = 0; index < Count; index++)
+            for (var index = 0; index < this.Count; index++)
             {
                 stringBuilder.Append(this[index].ToString(format, formatProvider));
-                if (index != Count - 1)
+                if (index != this.Count - 1)
                 {
                     stringBuilder.Append(formatProvider.GetTextInfo().ListSeparator);
                 }
@@ -1050,13 +1236,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Determines whether the specified <see cref="System.Object"/> is equal to this instance.
         /// </summary>
-        /// <param name="obj">The <see cref="System.Object"/> to compare with this instance.</param>
+        /// <param name="obj">
+        /// The <see cref="System.Object"/> to compare with this instance.
+        /// </param>
         /// <returns>
-        ///     <c>true</c> if the specified <see cref="System.Object"/> is equal to this instance; otherwise, <c>false</c>.
+        /// <c>true</c> if the specified <see cref="System.Object"/> is equal to this instance; otherwise, <c>false</c>.
         /// </returns>
         public override bool Equals(object obj)
         {
-            return Equals(obj as Vector);
+            return this.Equals(obj as Vector);
         }
 
         /// <summary>
@@ -1067,7 +1255,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public override int GetHashCode()
         {
-            var hashNum = Math.Min(Count, 20);
+            var hashNum = Math.Min(this.Count, 20);
             long hash = 0;
             for (var i = 0; i < hashNum; i++)
             {
@@ -1089,8 +1277,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public override string ToString()
         {
-            return ToString(null, null);
+            return this.ToString(null, null);
         }
+
         #endregion
     }
 }
\ No newline at end of file
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 1bac2a3e..4e20e60e 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -180,11 +180,6 @@
     <Compile Include="Statistics\MCMC\RejectionSampler.cs" />
     <Compile Include="Statistics\Statistics.cs" />
     <Compile Include="Statistics\MCMC\UnivariateSliceSampler.cs" />
-    <Compile Include="Threading\AggregateException.cs" />
-    <Compile Include="Threading\TaskOfT.cs" />
-    <Compile Include="Threading\Parallel.cs" />
-    <Compile Include="Threading\Task.cs" />
-    <Compile Include="Threading\ThreadQueue.cs" />
     <Compile Include="Trigonometry.cs" />
   </ItemGroup>
   <ItemGroup>
diff --git a/src/Numerics/Threading/AggregateException.cs b/src/Numerics/Threading/AggregateException.cs
deleted file mode 100644
index b376f19f..00000000
--- a/src/Numerics/Threading/AggregateException.cs
+++ /dev/null
@@ -1,69 +0,0 @@
-// <copyright file="AggregateException.cs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
-// Permission is hereby granted, free of charge, to any person
-// obtaining a copy of this software and associated documentation
-// files (the "Software"), to deal in the Software without
-// restriction, including without limitation the rights to use,
-// copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following
-// conditions:
-//
-// The above copyright notice and this permission notice shall be
-// included in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-// OTHER DEALINGS IN THE SOFTWARE.
-// </copyright>
-
-namespace MathNet.Numerics.Threading
-{
-    using System;
-    using System.Collections.Generic;
-    using System.Collections.ObjectModel;
-
-    /// <summary>
-    /// Represents multiple errors that occur during application execution.
-    /// </summary>
-#if !SILVERLIGHT
-    [Serializable]
-#endif
-    public class AggregateException : Exception
-    {
-        /// <summary>
-        /// List of the aggregated exceptions.
-        /// </summary>
-        private readonly IList<Exception> _exceptions = new List<Exception>();
-
-        /// <summary>
-        /// Initializes a new instance of the AggregateException class with a specified error message and references to the inner exceptions that are the cause of this exception.
-        /// </summary>
-        /// <param name="exceptions">The exceptions that are the cause of the current exception.</param>
-        public AggregateException(IEnumerable<Exception> exceptions)
-        {
-            foreach (var exception in exceptions)
-            {
-                _exceptions.Add(exception);
-            }
-        }
-
-        /// <summary>
-        /// Gets a read-only collection of the Exception instances that caused the current exception. 
-        /// </summary>
-        /// <value>A read-only collection of the Exception instances that caused the current exception</value>
-        public ReadOnlyCollection<Exception> InnerExceptions
-        {
-            get { return new ReadOnlyCollection<Exception>(_exceptions); }
-        }
-    }
-}
diff --git a/src/Numerics/Threading/Parallel.cs b/src/Numerics/Threading/Parallel.cs
deleted file mode 100644
index 7d87b68d..00000000
--- a/src/Numerics/Threading/Parallel.cs
+++ /dev/null
@@ -1,485 +0,0 @@
-// <copyright file="Parallel.cs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
-// Permission is hereby granted, free of charge, to any person
-// obtaining a copy of this software and associated documentation
-// files (the "Software"), to deal in the Software without
-// restriction, including without limitation the rights to use,
-// copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following
-// conditions:
-//
-// The above copyright notice and this permission notice shall be
-// included in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-// OTHER DEALINGS IN THE SOFTWARE.
-// </copyright>
-
-namespace MathNet.Numerics.Threading
-{
-    using System;
-    using System.Collections.Generic;
-    using Properties;
-
-    /// <summary>
-    /// Provides support for parallel loops.
-    /// </summary>
-    internal static class Parallel
-    {
-        /// <summary>
-        /// The amount to scale the foreach buffer after each iteration.
-        /// </summary>
-        private const int ScalingFactor = 2;
-
-        /// <summary>
-        /// The maximum size of the foreach buffer.
-        /// </summary>
-        private const int MaxBlockSize = 65536;
-
-        /// <summary>
-        /// The initial size of the for each buffer.
-        /// </summary>
-        private const int IntialBlockSize = 1024;
-
-        /// <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>
-        /// <exception cref="ArgumentNullException">The <paramref name="body"/> argument is null.</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)
-        {
-            if (body == null)
-            {
-                throw new ArgumentNullException("body");
-            }
-
-            // fast forward execution if it's only one or none items
-            var count = toExclusive - fromInclusive;
-            if (count <= 1)
-            {
-                if (count == 1)
-                {
-                    body(fromInclusive);
-                }
-
-                return;
-            }
-
-            // fast forward execution in case parallelization is disabled
-            if (Control.DisableParallelization
-                || ThreadQueue.ThreadCount <= 1
-                || ThreadQueue.IsInWorkerThread)
-            {
-                for (int i = fromInclusive; i < toExclusive; i++)
-                {
-                    body(i);
-                }
-
-                return;
-            }
-
-            var actions = new Action[ThreadQueue.ThreadCount];
-            var size = count / actions.Length;
-
-            // partition the jobs into separate sets for each but the last worked thread
-            for (var i = 0; i < actions.Length - 1; i++)
-            {
-                var start = fromInclusive + (i * size);
-                var stop = fromInclusive + ((i + 1) * size);
-
-                actions[i] =
-                    () =>
-                    {
-                        for (int j = start; j < stop; j++)
-                        {
-                            body(j);
-                        }
-                    };
-            }
-
-            // add another set for last worker thread
-            actions[actions.Length - 1] =
-                () =>
-                {
-                    for (int i = fromInclusive + ((actions.Length - 1) * size); i < toExclusive; i++)
-                    {
-                        body(i);
-                    }
-                };
-
-            Invoke(actions);
-        }
-
-        /// <summary>
-        /// Executes a for loop in which iterations may run in parallel.
-        /// </summary>
-        /// <typeparam name="T">The type of the thread-local data.</typeparam>
-        /// <param name="fromInclusive">The start index, inclusive.</param>
-        /// <param name="toExclusive">The end index, exclusive.</param>
-        /// <param name="localInit">The function delegate that returns the initial state of the local data for each thread.</param>
-        /// <param name="body">The delegate that is invoked once per iteration.</param>
-        /// <param name="localFinally">The delegate that performs a final action on the local state of each thread.</param>
-        public static void For<T>(int fromInclusive, int toExclusive, Func<T> localInit, Func<int, T, T> body, Action<T> localFinally)
-        {
-            var count = toExclusive - fromInclusive;
-            var tasks = new Task<T>[ThreadQueue.ThreadCount];
-            var size = count / tasks.Length;
-
-            // fast forward execution if it's only one or none items            
-            if (count <= 1)
-            {
-                if (count == 1)
-                {
-                    localFinally(body(fromInclusive, localInit()));
-                }
-
-                return;
-            }
-
-            // fast forward execution in case parallelization is disabled
-            if (Control.DisableParallelization
-                || ThreadQueue.ThreadCount <= 1
-                || ThreadQueue.IsInWorkerThread)
-            {
-                var localresult = localInit();
-                for (var i = fromInclusive; i < toExclusive; i++)
-                {
-                    localresult = body(i, localresult);
-                }
-
-                localFinally(localresult);
-                return;
-            }
-
-            // partition the jobs into separate sets for each but the last worked thread
-            for (var i = 0; i < tasks.Length - 1; i++)
-            {
-                var start = fromInclusive + (i * size);
-                var stop = fromInclusive + ((i + 1) * size);
-                tasks[i] = new Task<T>(
-                    localData =>
-                    {
-                        var localresult = (T)localData;
-                        for (var j = start; j < stop; j++)
-                        {
-                            localresult = body(j, localresult);
-                        }
-
-                        return localresult;
-                    },
-                    localInit());
-                ThreadQueue.Enqueue(tasks[i]);
-            }
-
-            // add another set for last worker thread
-            tasks[tasks.Length - 1] = new Task<T>(
-                    localData =>
-                    {
-                        var localresult = (T)localData;
-                        for (var i = fromInclusive + ((tasks.Length - 1) * size); i < toExclusive; i++)
-                        {
-                            localresult = body(i, localresult);
-                        }
-
-                        return localresult;
-                    }, 
-                    localInit());
-
-            ThreadQueue.Enqueue(tasks[tasks.Length - 1]);
-            if (tasks.Length <= 0)
-            {
-                return;
-            }
-
-            WaitForTasksToComplete(tasks);
-
-            foreach (var t in tasks)
-            {
-                localFinally(t.Result);
-            }
-
-            CollectExceptions(tasks);
-        }
-
-        /// <summary>
-        /// Executes a for each operation on an IEnumerable{T} in which iterations may run in parallel.
-        /// </summary>
-        /// <typeparam name="T">The type of the data in the source.</typeparam>
-        /// <param name="source">An enumerable data source.</param>
-        /// <param name="body">The delegate that is invoked once per iteration.</param>
-        public static void ForEach<T>(IEnumerable<T> source, Action<T> body)
-        {
-            if (body == null)
-            {
-                throw new ArgumentNullException("body");
-            }
-
-            // fast forward execution in case parallelization is disabled
-            if (Control.DisableParallelization
-                || ThreadQueue.ThreadCount <= 1
-                || ThreadQueue.IsInWorkerThread)
-            {
-                foreach (var item in source)
-                {
-                    body(item);
-                }
-
-                return;
-            }
-
-            // source is a IList, call For instead.
-            if (source is IList<T>)
-            {
-                var list = (IList<T>)source;
-                For(0, list.Count, i => body(list[i]));
-                return;
-            }
-
-            var maxBlockSize = IntialBlockSize;
-            var tasks = new List<Task>();
-
-            var enumerator = source.GetEnumerator(); 
-            while (enumerator.MoveNext())
-            {
-                var pos = 0;
-                var list = new T[maxBlockSize];
-                list[pos++] = enumerator.Current;
-
-                var count = 1;
-                while (count < maxBlockSize && enumerator.MoveNext())
-                {
-                    list[pos++] = enumerator.Current;
-                    count++;
-                }
-
-                var task = new Task(
-                    () =>
-                    {
-                        for (var i = 0; i < pos; i++)
-                        {
-                            body(list[i]);
-                        }
-                    });
-
-                ThreadQueue.Enqueue(task);
-                tasks.Add(task);
-                maxBlockSize = Math.Min(MaxBlockSize, maxBlockSize * ScalingFactor);
-            }
-
-            if (tasks.Count > 0)
-            {
-                WaitForTasksToComplete(tasks.ToArray());
-                CollectExceptions(tasks);
-            }
-        }
-
-        /// <summary>
-        /// Executes a for each operation on an IEnumerable{TSource in which iterations may run in parallel.
-        /// </summary>
-        /// <typeparam name="TSource">The type of the data in the source.</typeparam>
-        /// <typeparam name="TLocal">The type of the thread-local data.</typeparam>
-        /// <param name="source">An enumerable data source.</param>
-        /// <param name="localInit">The function delegate that returns the initial state of the local data for each thread.</param>
-        /// <param name="body">The delegate that is invoked once per iteration.</param>
-        /// <param name="localFinally">The delegate that performs a final action on the local state of each thread.</param>
-        public static void ForEach<TSource, TLocal>(IEnumerable<TSource> source, Func<TLocal> localInit, Func<TSource, TLocal, TLocal> body, Action<TLocal> localFinally)
-        {
-            if (body == null)
-            {
-                throw new ArgumentNullException("body");
-            }
-
-            // fast forward execution in case parallelization is disabled
-            if (Control.DisableParallelization
-                || ThreadQueue.ThreadCount <= 1
-                || ThreadQueue.IsInWorkerThread)
-            {
-                var localResult = localInit();
-                foreach (var item in source)
-                {
-                    localResult = body(item, localResult);
-                }
-
-                localFinally(localResult);
-                return;
-            }
-
-            // source is a IList, call For instead.
-            if (source is IList<TSource>)
-            {
-                var list = (IList<TSource>)source;
-                For(0, list.Count, localInit, (i, local) => body(list[i], local), localFinally);
-                return;
-            }
-
-            var maxBlockSize = IntialBlockSize;
-            var tasks = new List<Task<TLocal>>();
-
-            var enumerator = source.GetEnumerator();
-            while (enumerator.MoveNext())
-            {
-                var pos = 0;
-                var list = new TSource[maxBlockSize];
-                list[pos++] = enumerator.Current;
-
-                var count = 1;
-                while (count < maxBlockSize && enumerator.MoveNext())
-                {
-                    list[pos++] = enumerator.Current;
-                    count++;
-                }
-
-                var task = new Task<TLocal>(
-                    localData =>
-                    {
-                        var localresult = localData;
-                        for (var i = 0; i < pos; i++)
-                        {
-                            localresult = body(list[i], (TLocal)localresult);
-                        }
-
-                        return (TLocal)localresult;
-                    }, 
-                    localInit());
-
-                ThreadQueue.Enqueue(task);
-                tasks.Add(task);
-                maxBlockSize = Math.Min(MaxBlockSize, maxBlockSize * ScalingFactor);
-            }
-
-            if (tasks.Count <= 0)
-            {
-                return;
-            }
-            
-            var taskArray = tasks.ToArray();
-            WaitForTasksToComplete(taskArray);
-
-            for (var i = 0; i < taskArray.Length; i++)
-            {
-                localFinally(tasks[i].Result);
-            }
-
-            CollectExceptions(taskArray);
-        }
-
-        /// <summary>
-        /// Executes each of the provided actions inside a discrete, asynchronous task. 
-        /// </summary>
-        /// <param name="actions">An array of actions to execute.</param>
-        /// <exception cref="ArgumentNullException">The <paramref name="actions"/> argument is null.</exception>
-        /// <exception cref="ArgumentException">The actions array contains a null element.</exception>
-        /// <exception cref="AggregateException">An action threw an exception.</exception>
-        public static void Run(params Action[] actions)
-        {
-            if (actions == null)
-            {
-                throw new ArgumentNullException("actions");
-            }
-
-            // fast forward execution if it's only one or none items
-            if (actions.Length <= 1)
-            {
-                if (actions.Length == 1)
-                {
-                    actions[0]();
-                }
-
-                return;
-            }
-
-            // fast forward execution in case parallelization is disabled
-            if (Control.DisableParallelization
-                || ThreadQueue.ThreadCount <= 1
-                || ThreadQueue.IsInWorkerThread)
-            {
-                for (var i = 0; i < actions.Length; i++)
-                {
-                    actions[i]();
-                }
-
-                return;
-            }
-
-            Invoke(actions);
-        }
-
-        /// <summary>
-        /// Executes each of the provided actions inside a discrete, asynchronous task. 
-        /// </summary>
-        /// <param name="actions">An array of actions to execute.</param>
-        /// <exception cref="ArgumentException">The actions array contains a null element.</exception>
-        /// <exception cref="AggregateException">An action threw an exception.</exception>
-        private static void Invoke(params Action[] actions)
-        {
-            // create a job for each action
-            var tasks = new Task[actions.Length];
-            for (var i = 0; i < tasks.Length; i++)
-            {
-                Action action = actions[i];
-                if (action == null)
-                {
-                    throw new ArgumentException(String.Format(Resources.ArgumentItemNull, "actions"), "actions");
-                }
-
-                tasks[i] = new Task(action);
-            }
-
-            // run the jobs
-            ThreadQueue.Enqueue(tasks);
-
-            WaitForTasksToComplete(tasks);
-
-            CollectExceptions(tasks);
-        }
-
-        /// <summary>
-        /// Waits for tasks to complete.
-        /// </summary>
-        /// <param name="tasks">The tasks.</param>
-        private static void WaitForTasksToComplete(Task[] tasks)
-        {
-            foreach (var task in tasks)
-            {
-                task.Wait();
-            }
-        }
-
-        /// <summary>
-        /// Collects the exceptions and dispose tasks.
-        /// </summary>
-        /// <param name="tasks">The tasks.</param>
-        private static void CollectExceptions(IEnumerable<Task> tasks)
-        {
-            // collect all thrown exceptions and dispose the jobs
-            var exceptions = new List<Exception>();
-            foreach (var task in tasks)
-            {
-                if (task.IsFaulted)
-                {
-                    exceptions.Add(task.Exception);
-                }
-            }
-
-            // throw the aggregated exceptions, if any
-            if (exceptions.Count > 0)
-            {
-                throw new AggregateException(exceptions);
-            }
-        }
-    }
-}
diff --git a/src/Numerics/Threading/Task.cs b/src/Numerics/Threading/Task.cs
deleted file mode 100644
index 97ffd43a..00000000
--- a/src/Numerics/Threading/Task.cs
+++ /dev/null
@@ -1,127 +0,0 @@
-// <copyright file="Task.cs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
-// Permission is hereby granted, free of charge, to any person
-// obtaining a copy of this software and associated documentation
-// files (the "Software"), to deal in the Software without
-// restriction, including without limitation the rights to use,
-// copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following
-// conditions:
-//
-// The above copyright notice and this permission notice shall be
-// included in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-// OTHER DEALINGS IN THE SOFTWARE.
-// </copyright>
-
-namespace MathNet.Numerics.Threading
-{
-    using System;
-    using System.Threading;
-
-    /// <summary>
-    /// Internal Parallel Task Handle.
-    /// </summary>
-    internal class Task
-    {
-        /// <summary>
-        /// Delegate to the task's action.
-        /// </summary>
-        private readonly Action _body;
-
-        /// <summary>
-        /// Initializes a new instance of the Task class.
-        /// </summary>
-        /// <param name="body">Delegate to the task's action.</param>
-        public Task(Action body)
-        {
-            if (body == null)
-            {
-                throw new ArgumentNullException("body");
-            }
-
-            _body = body;
-        }
-
-        /// <summary>
-        /// Initializes a new instance of the <see cref="Task"/> class.
-        /// </summary>
-        protected Task()
-        {
-        }
-
-        /// <summary>
-        /// Gets a value indicating whether the task completed due to an unhandled exception.
-        /// </summary>
-        /// <value>
-        ///     <c>true</c> if this task completed due to an unhandled exception; otherwise, <c>false</c>.
-        /// </value>
-        public bool IsFaulted
-        {
-            get { return Exception != null; }
-        }
-
-        /// <summary>
-        /// Gets a value indicating whether this task has completed.
-        /// </summary>
-        /// <value>
-        ///     <c>true</c> if this task has completed; otherwise, <c>false</c>.
-        /// </value>
-        public bool IsCompleted
-        {
-            get; private set;
-        }
-
-        /// <summary>
-        /// Gets or sets the exception thrown by the task, if any.
-        /// </summary>
-        public Exception Exception { get; set; }
-
-        /// <summary>
-        /// Run the task.
-        /// </summary>
-        public void Compute()
-        {
-            try
-            {
-                DoCompute();
-                IsCompleted = true;
-            }
-            catch (Exception e)
-            {
-                Exception = e;
-            }
-        }
-
-        /// <summary>
-        /// Runs the actual task.
-        /// </summary>
-        protected virtual void DoCompute()
-        {
-            _body();
-        }
-
-        /// <summary>
-        /// Waits for the task to complete execution.
-        /// </summary>
-        public void Wait()
-        {
-            while (!IsCompleted && !IsFaulted)
-            {
-                Thread.Sleep(0);
-            }
-        }
-    }
-}
diff --git a/src/Numerics/Threading/TaskOfT.cs b/src/Numerics/Threading/TaskOfT.cs
deleted file mode 100644
index 9638e874..00000000
--- a/src/Numerics/Threading/TaskOfT.cs
+++ /dev/null
@@ -1,79 +0,0 @@
-// <copyright file="TaskOfT.cs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
-// Permission is hereby granted, free of charge, to any person
-// obtaining a copy of this software and associated documentation
-// files (the "Software"), to deal in the Software without
-// restriction, including without limitation the rights to use,
-// copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following
-// conditions:
-//
-// The above copyright notice and this permission notice shall be
-// included in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-// OTHER DEALINGS IN THE SOFTWARE.
-// </copyright>
-
-namespace MathNet.Numerics.Threading
-{
-    using System;
-
-    /// <summary>
-    /// Internal Generic Parallel Task Handle.
-    /// </summary>
-    /// <typeparam name="TResult">The type of the result.</typeparam>
-    internal class Task<TResult> : Task
-    {
-        /// <summary>
-        /// Delegate to the task's action.
-        /// </summary>
-        private readonly Func<object, TResult> _body;
-
-        /// <summary>
-        /// Variable used to hold state information between iterations.
-        /// </summary>
-        private readonly object _state;
-        
-        /// <summary>
-        /// Gets the result of the task.
-        /// </summary>
-        /// <value>The result of the task.</value>
-        public TResult Result { get; private set; }
-
-        /// <summary>
-        /// Initializes a new instance of the Task class.
-        /// </summary>
-        /// <param name="body">Delegate to the task's action.</param>
-        /// <param name="state">An object representing data to be used by the action.</param>
-        public Task(Func<object, TResult> body, object state)
-        {
-            if (body == null)
-            {
-                throw new ArgumentNullException("body");
-            }
-
-            _state = state;
-            _body = body;
-        }
-
-        /// <summary>
-        /// Runs the actual task.
-        /// </summary>
-        protected override void DoCompute()
-        {
-            Result = _body(_state);
-        }
-    }
-}
diff --git a/src/Numerics/Threading/ThreadQueue.cs b/src/Numerics/Threading/ThreadQueue.cs
deleted file mode 100644
index 41320536..00000000
--- a/src/Numerics/Threading/ThreadQueue.cs
+++ /dev/null
@@ -1,262 +0,0 @@
-// <copyright file="ThreadQueue.cs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
-// Permission is hereby granted, free of charge, to any person
-// obtaining a copy of this software and associated documentation
-// files (the "Software"), to deal in the Software without
-// restriction, including without limitation the rights to use,
-// copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following
-// conditions:
-//
-// The above copyright notice and this permission notice shall be
-// included in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-// OTHER DEALINGS IN THE SOFTWARE.
-// </copyright>
-
-namespace MathNet.Numerics.Threading
-{
-    using System;
-    using System.Collections.Generic;
-    using System.Threading;
-
-    /// <summary>
-    /// Internal Parallel Thread Queue.
-    /// </summary>
-    internal static class ThreadQueue
-    {
-        /// <summary>
-        /// Sync Object for the thread queue state.
-        /// </summary>
-        private static readonly object _stateSync = new object();
-
-        /// <summary>
-        /// Sync Object for queue access (to be sure it's used by us only).
-        /// </summary>
-        private static readonly object _queueSync = new object();
-
-        /// <summary>
-        /// Queue holding the pending jobs.
-        /// </summary>
-        private static readonly Queue<Task> _queue = new Queue<Task>();
-
-        /// <summary>
-        /// Running flag, used to signal worker threads to stop cleanly.
-        /// </summary>
-        private static bool _running = true;
-
-        /// <summary>
-        /// Worker threads
-        /// </summary>
-        private static Thread[] _threads;
-
-        /// <summary>
-        /// Gets the number of worker threads.
-        /// </summary>
-        internal static int ThreadCount { get; private set; }
-
-        /// <summary>
-        /// Indicating whether the current thread is a parallelized worker thread.
-        /// </summary>
-        [ThreadStatic]
-        private static bool _isInWorkerThread;
-
-        /// <summary>
-        /// Initializes static members of the ThreadQueue class.
-        /// </summary>
-        static ThreadQueue()
-        {
-            #if SILVERLIGHT
-            Start(1);
-            #else
-            Start(Environment.ProcessorCount);
-            #endif
-        }
-
-        /// <summary>
-        /// Gets a value indicating whether the current thread is a parallelized worker thread.
-        /// </summary>
-        public static bool IsInWorkerThread
-        {
-            get { return _isInWorkerThread; }
-        }
-
-        /// <summary>
-        /// Add a job to the queue.
-        /// </summary>
-        /// <param name="task">The job to run.</param>
-        public static void Enqueue(Task task)
-        {
-            if (!_running)
-            {
-                Start();
-            }
-
-            lock (_queueSync)
-            {
-                _queue.Enqueue(task);
-                Monitor.Pulse(_queueSync);
-            }
-        }
-
-        /// <summary>
-        /// Add a set of jobs to the queue.
-        /// </summary>
-        /// <param name="tasks">The jobs to run.</param>
-        public static void Enqueue(IList<Task> tasks)
-        {
-            if (!_running)
-            {
-                Start();
-            }
-
-            lock (_queueSync)
-            {
-                foreach (var task in tasks)
-                {
-                    _queue.Enqueue(task);
-                }
-
-                Monitor.PulseAll(_queueSync);
-            }
-        }
-
-        /// <summary>
-        /// Worker Thread Program
-        /// </summary>
-        private static void WorkerThreadStart()
-        {
-            _isInWorkerThread = true;
-
-            while (_running)
-            {
-                // Get the job...
-                Task task = null;
-
-                lock (_queueSync)
-                {
-                    // Check whether we should shut down
-                    if (!_running)
-                    {
-                        break;
-                    }
-
-                    if (_queue.Count > 0)
-                    {
-                        task = _queue.Dequeue();
-                    }
-                    else
-                    {
-                        Monitor.Wait(_queueSync);
-                    }
-                }
-
-                if (task == null)
-                {
-                    continue;
-                }
-
-                // ...and run it
-                task.Compute();
-            }
-
-            _isInWorkerThread = false;
-        }
-
-        /// <summary>
-        /// Start or restart the queue with the specified number of worker threads.
-        /// </summary>
-        /// <param name="numberOfThreads">Number of worker threads.</param>
-        public static void Start(int numberOfThreads)
-        {
-            lock (_stateSync)
-            {
-                // instead of throwing an out of range exception, simply normalize
-                numberOfThreads = Math.Max(1, Math.Min(1024, numberOfThreads));
-
-                if (_threads != null)
-                {
-                    if (_threads.Length == numberOfThreads)
-                    {
-                        return;
-                    }
-
-                    Shutdown();
-                }
-
-                ThreadCount = numberOfThreads;
-                Start();
-            }
-         }
-
-        /// <summary>
-        /// Start the thread queue, if it is not already running.
-        /// </summary>
-        public static void Start()
-        {
-            lock (_stateSync)
-            {
-                if (_threads != null)
-                {
-                    return;
-                }
-
-                _running = true; 
-                
-                _threads = new Thread[ThreadCount];
-
-                for (var i = 0; i < _threads.Length; i++)
-                {
-                    _threads[i] = new Thread(WorkerThreadStart)
-                                  {
-                                      IsBackground = true
-                                  };
-
-                    _threads[i].Start();
-                }
-            }
-        }
-
-        /// <summary>
-        /// Stop the thread queue, if it is running.
-        /// </summary>
-        public static void Shutdown()
-        {
-            // try to stop the worker threads cleanly
-            lock (_stateSync)
-            {
-                if (_threads == null)
-                {
-                    return;
-                }
-                
-                _running = false;
-
-                lock (_queueSync)
-                {
-                    Monitor.PulseAll(_queueSync);
-                }
-
-                // wait until all threads have stopped
-                foreach (var thread in _threads)
-                {
-                    thread.Join();
-                }
-                
-                _threads = null;
-            }
-        }
-    }
-}
diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj
index eacf7d5a..e5b3abf5 100644
--- a/src/Silverlight/Silverlight.csproj
+++ b/src/Silverlight/Silverlight.csproj
@@ -338,21 +338,6 @@
     <Compile Include="..\Numerics\Statistics\Statistics.cs">
       <Link>Statistics\Statistics.cs</Link>
     </Compile>
-    <Compile Include="..\Numerics\Threading\AggregateException.cs">
-      <Link>Threading\AggregateException.cs</Link>
-    </Compile>
-    <Compile Include="..\Numerics\Threading\Parallel.cs">
-      <Link>Threading\Parallel.cs</Link>
-    </Compile>
-    <Compile Include="..\Numerics\Threading\Task.cs">
-      <Link>Threading\Task.cs</Link>
-    </Compile>
-    <Compile Include="..\Numerics\Threading\TaskOfT.cs">
-      <Link>Threading\TaskOfT.cs</Link>
-    </Compile>
-    <Compile Include="..\Numerics\Threading\ThreadQueue.cs">
-      <Link>Threading\ThreadQueue.cs</Link>
-    </Compile>
     <Compile Include="..\Numerics\Trigonometry.cs">
       <Link>Trigonometry.cs</Link>
     </Compile>
diff --git a/src/UnitTests/ThreadingTests/ParallelForEachTests.cs b/src/UnitTests/ThreadingTests/ParallelForEachTests.cs
deleted file mode 100644
index e8a48f19..00000000
--- a/src/UnitTests/ThreadingTests/ParallelForEachTests.cs
+++ /dev/null
@@ -1,602 +0,0 @@
-// <copyright file="ParallelTest.cs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
-// Permission is hereby granted, free of charge, to any person
-// obtaining a copy of this software and associated documentation
-// files (the "Software"), to deal in the Software without
-// restriction, including without limitation the rights to use,
-// copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following
-// conditions:
-//
-// The above copyright notice and this permission notice shall be
-// included in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-// OTHER DEALINGS IN THE SOFTWARE.
-// </copyright>
-
-namespace MathNet.Numerics.UnitTests.ThreadingTests
-{
-    using System;
-    using System.Collections.Generic;
-    using System.Threading;
-    using MbUnit.Framework;
-    using Threading;
-
-    [TestFixture]
-    public class ParallelForEachTests
-    {
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceMTAOnePerCore(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceSTAOnePerCore(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceMTATwoPerCore(int count)
-        {
-            // ensure Two-Per-Core
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceSTATwoPerCore(int count)
-        {
-            // ensure Two-Per-Core
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceMTAOnePerCore(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceSTAOnePerCore(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceMTATwoPerCore(int count)
-        {
-            // ensure Two-Per-Core
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceSTATwoPerCore(int count)
-        {
-            // ensure Two-Per-Core
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var items = new double[count];
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-            Parallel.ForEach(pairs, pair => items[pair.Key] = pair.Value);
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        public void DoesNotGetConfusedByMultipleStartShutdown()
-        {
-            ThreadQueue.Shutdown();
-            ThreadQueue.Shutdown();
-
-            ThreadQueue.Start(2);
-            Assert.AreEqual(2, ThreadQueue.ThreadCount);
-
-            Control.NumberOfParallelWorkerThreads = 2;
-            Assert.AreEqual(2, ThreadQueue.ThreadCount);
-
-            ThreadQueue.Start(4);
-            Assert.AreEqual(4, ThreadQueue.ThreadCount);
-            Assert.AreEqual(4, Control.NumberOfParallelWorkerThreads);
-
-            ThreadQueue.Shutdown();
-            ThreadQueue.Start();
-            Assert.AreEqual(4, ThreadQueue.ThreadCount);
-
-            ThreadQueue.Start(2);
-            Assert.AreEqual(2, ThreadQueue.ThreadCount);
-
-            var items = new double[50];
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < items.Length; i++)
-            {
-                items[i] = i;
-                pairs.Add(new KeyValuePair<int, double>(i, i));
-            }
-
-            Parallel.ForEach(pairs,
-                pair =>
-                {
-                    items[pair.Key] = pair.Value + 1000;
-                }
-                );
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1000 + i, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        public void DoesDetectAndResolveRecursiveParallelization()
-        {
-            var countSharedBetweenClosures = 0;
-            var values = new int[10];
-
-            Assert.DoesNotThrow(
-                () =>
-                Parallel.ForEach(values,
-                    j => Interlocked.Increment(ref countSharedBetweenClosures)));
-
-            Assert.AreEqual(10, countSharedBetweenClosures);
-            countSharedBetweenClosures = 0;
-
-            Parallel.ForEach(values,
-                i =>
-                Parallel.ForEach(values,
-                    j => Interlocked.Increment(ref countSharedBetweenClosures)));
-
-            Assert.AreEqual(100, countSharedBetweenClosures);
-            countSharedBetweenClosures = 0;
-
-            Parallel.ForEach(values,
-                i =>
-                Parallel.ForEach(values,
-                    j =>
-                    Parallel.ForEach(values,
-                        k => Interlocked.Increment(ref countSharedBetweenClosures))));
-
-            Assert.AreEqual(1000, countSharedBetweenClosures);
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceMTAOnePerCoreWithIntialAndFinally(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceSTAOnePerCoreWithIntialAndFinally(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceMTATwoPerCoreWithIntialAndFinally(int count)
-        {
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInListOnceSTATwoPerCoreWithIntialAndFinally(int count)
-        {
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var pairs = new List<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceMTAOnePerCoreWithIntialAndFinally(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceSTAOnePerCoreWithIntialAndFinally(int count)
-        {
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceMTATwoPerCoreWithIntialAndFinally(int count)
-        {
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1001)]
-        public void ParallelForEachInvokesEveryItemInSetOnceSTATwoPerCoreWithIntialAndFinally(int count)
-        {
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var pairs = new HashSet<KeyValuePair<int, double>>();
-            for (var i = 0; i < count; i++)
-            {
-                pairs.Add(new KeyValuePair<int, double>(i, 2));
-            }
-
-            var sum = 0.0;
-            var sync = new object();
-
-            Parallel.ForEach(
-                pairs, () => 0.0, (pair, localData) => localData += pair.Value,
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, Timeout(15)]
-        public void DoesDetectAndResolveRecursiveParallelizationWithIntialAndFinally()
-        {
-            var countSharedBetweenClosures = 0;
-            var values = new int[10];
-
-            Assert.DoesNotThrow(
-                () =>
-                Parallel.ForEach(values,
-                    j => Interlocked.Increment(ref countSharedBetweenClosures)));
-
-            Assert.AreEqual(10, countSharedBetweenClosures);
-            countSharedBetweenClosures = 0;
-
-            Parallel.ForEach(values,
-                i =>
-                Parallel.ForEach(values,
-                    j => Interlocked.Increment(ref countSharedBetweenClosures)));
-
-            Assert.AreEqual(100, countSharedBetweenClosures);
-            countSharedBetweenClosures = 0;
-
-            Parallel.ForEach(values,
-                i =>
-                Parallel.ForEach(values,
-                    j =>
-                    Parallel.ForEach(values, () => 0.0,
-                        (iterator, localData) => Interlocked.Increment(ref countSharedBetweenClosures),
-                        localResult => { return; }
-                        )));
-
-            Assert.AreEqual(1000, countSharedBetweenClosures);
-        }
-    }
-}
\ No newline at end of file
diff --git a/src/UnitTests/ThreadingTests/ParallelForTests.cs b/src/UnitTests/ThreadingTests/ParallelForTests.cs
deleted file mode 100644
index ea0c344c..00000000
--- a/src/UnitTests/ThreadingTests/ParallelForTests.cs
+++ /dev/null
@@ -1,364 +0,0 @@
-// <copyright file="ParallelTest.cs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://mathnet.opensourcedotnet.info
-//
-// Copyright (c) 2009 Math.NET
-//
-// Permission is hereby granted, free of charge, to any person
-// obtaining a copy of this software and associated documentation
-// files (the "Software"), to deal in the Software without
-// restriction, including without limitation the rights to use,
-// copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following
-// conditions:
-//
-// The above copyright notice and this permission notice shall be
-// included in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-// OTHER DEALINGS IN THE SOFTWARE.
-// </copyright>
-
-namespace MathNet.Numerics.UnitTests.ThreadingTests
-{
-    using System;
-    using System.Threading;
-    using MbUnit.Framework;
-    using Threading;
-
-    [TestFixture]
-    public class ParallelForTests
-    {
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceMTAOnePerCore(int count)
-        {
-            var items = new int[count];
-
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            Parallel.For(0, count, i => items[i]++);
-            Parallel.For(0, count, i => items[i] += 1000);
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1001, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceSTAOnePerCore(int count)
-        {
-            var items = new int[count];
-
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            Parallel.For(0, count, i => items[i]++);
-            Parallel.For(0, count, i => items[i] += 1000);
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1001, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceMTATwoPerCore(int count)
-        {
-            var items = new int[count];
-
-            // ensure Two-Per-Core
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            Parallel.For(0, count, i => items[i]++);
-            Parallel.For(0, count, i => items[i] += 1000);
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1001, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceSTATwoPerCore(int count)
-        {
-            var items = new int[count];
-
-            // ensure Two-Per-Core
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            Parallel.For(0, count, i => items[i]++);
-            Parallel.For(0, count, i => items[i] += 1000);
-
-            for (int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1001, items[i], i.ToString());
-            }
-        }
-
-        [Test, MultipleAsserts, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        public void DoesNotGetConfusedByMultipleStartShutdown()
-        {
-            ThreadQueue.Shutdown();
-            ThreadQueue.Shutdown();
-
-            ThreadQueue.Start(2);
-            Assert.AreEqual(2, ThreadQueue.ThreadCount);
-
-            Control.NumberOfParallelWorkerThreads = 2;
-            Assert.AreEqual(2, ThreadQueue.ThreadCount);
-
-            ThreadQueue.Start(4);
-            Assert.AreEqual(4, ThreadQueue.ThreadCount);
-            Assert.AreEqual(4, Control.NumberOfParallelWorkerThreads);
-
-            ThreadQueue.Shutdown();
-            ThreadQueue.Start();
-            Assert.AreEqual(4, ThreadQueue.ThreadCount);
-
-            ThreadQueue.Start(2);
-            Assert.AreEqual(2, ThreadQueue.ThreadCount);
-
-            var items = new int[50];
-            Parallel.For(0, items.Length, i => items[i]++);
-
-            Parallel.For(0, items.Length, i => items[i] += 1000);
-
-            ThreadQueue.Shutdown();
-
-            for(int i = 0; i < items.Length; i++)
-            {
-                Assert.AreEqual(1001, items[i], i.ToString());
-            }
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        public void DoesDetectAndResolveRecursiveParallelization()
-        {
-            int countSharedBetweenClosures = 0;
-
-            Assert.DoesNotThrow(
-                () =>
-                Parallel.For(
-                    0,
-                    10,
-                    j => Interlocked.Increment(ref countSharedBetweenClosures)));
-
-            Assert.AreEqual(10, countSharedBetweenClosures);
-            countSharedBetweenClosures = 0;
-
-            Parallel.For(
-                0,
-                10,
-                i =>
-                Parallel.For(
-                    0,
-                    10,
-                    j => Interlocked.Increment(ref countSharedBetweenClosures)));
-
-            Assert.AreEqual(100, countSharedBetweenClosures);
-            countSharedBetweenClosures = 0;
-
-            Parallel.For(
-                0,
-                10,
-                i =>
-                Parallel.For(
-                    0,
-                    10,
-                    j =>
-                    Parallel.For(
-                        0,
-                        10,
-                        k => Interlocked.Increment(ref countSharedBetweenClosures))));
-
-            Assert.AreEqual(1000, countSharedBetweenClosures);
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceMTAOnePerCoreWithIntialAndFinally(int count)
-        {
-            var items = new double[count];
-
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var sum = 0.0;
-            for(var i = 0; i < items.Length; i++)
-            {
-                items[i] = 2;
-            }
-            var sync = new object();
-            Parallel.For(
-                0,
-                count,
-                () => 0.0,
-                (i, localData) =>localData += items[i],
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceSTAOnePerCoreWithIntialAndFinally(int count)
-        {
-            var items = new double[count];
-
-            // ensure One-Per-Core
-            ThreadQueue.Start(Environment.ProcessorCount);
-
-            var sum = 0.0;
-            for (var i = 0; i < items.Length; i++)
-            {
-                items[i] = 2;
-            }
-            var sync = new object();
-            Parallel.For(
-                0,
-                count,
-                () => 0.0,
-                (i, localData) => localData += items[i],
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceMTATwoPerCoreWithIntialAndFinally(int count)
-        {
-            var items = new double[count];
-
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var sum = 0.0;
-            for (var i = 0; i < items.Length; i++)
-            {
-                items[i] = 2;
-            }
-            var sync = new object();
-            Parallel.For(
-                0,
-                count,
-                () => 0.0,
-                (i, localData) => localData += items[i],
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.STA), Timeout(15)]
-        [Column(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 101)]
-        public void ParallelForInvokesEveryItemOnceSTATwoPerCoreWithIntialAndFinally(int count)
-        {
-            var items = new double[count];
-
-            ThreadQueue.Start(2 * Environment.ProcessorCount);
-
-            var sum = 0.0;
-            for (var i = 0; i < items.Length; i++)
-            {
-                items[i] = 2;
-            }
-            var sync = new object();
-            Parallel.For(
-                0,
-                count,
-                () => 0.0,
-                (i, localData) => localData += items[i],
-                localResult =>
-                {
-                    lock (sync)
-                    {
-                        sum += localResult;
-                    }
-                }
-                );
-            Assert.AreEqual(count * 2, sum);
-        }
-
-        [Test, ApartmentState(ApartmentState.MTA), Timeout(15)]
-        public void DoesDetectAndResolveRecursiveParallelizationWithIntialAndFinally()
-       {
-           int countSharedBetweenClosures = 0;
-
-           Assert.DoesNotThrow(
-               () =>
-               Parallel.For(
-                   0,
-                   10,
-                   () => 0.0,
-                   (j, localData) => Interlocked.Increment(ref countSharedBetweenClosures),
-                   localResult => { return; }
-                   ));
-
-           Assert.AreEqual(10, countSharedBetweenClosures);
-           countSharedBetweenClosures = 0;
-
-           Parallel.For(
-               0,
-               10,
-               i =>
-               Parallel.For(
-                   0,
-                   10,
-                   () => 0.0,
-                   (j, localData) => Interlocked.Increment(ref countSharedBetweenClosures),
-                   localResult => { return; }
-                   ));
-
-           Assert.AreEqual(100, countSharedBetweenClosures);
-           countSharedBetweenClosures = 0;
-
-           Parallel.For(
-               0,
-               10,
-               i =>
-               Parallel.For(
-                   0,
-                   10,
-                   j =>
-                   Parallel.For(
-                       0,
-                       10,
-                       () => 0.0,
-                       (k, localData) => Interlocked.Increment(ref countSharedBetweenClosures),
-                       localResult => { return; }
-                       )));
-
-           Assert.AreEqual(1000, countSharedBetweenClosures);
-       }
-    }
-}
\ No newline at end of file
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index d2d15ea6..bcb623ce 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -153,8 +153,6 @@
     <Compile Include="StatisticsTests\StatisticsTests.cs" />
     <Compile Include="StatisticsTests\StatTestData.cs" />
     <Compile Include="StatisticsTests\MCMCTests\UnivariateSliceSamplerTests.cs" />
-    <Compile Include="ThreadingTests\ParallelForEachTests.cs" />
-    <Compile Include="ThreadingTests\ParallelForTests.cs" />
     <Compile Include="TrigonometryTest.cs" />
   </ItemGroup>
   <ItemGroup>