linear algebra: parallel dot product

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs

@@ -118,12 +118,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            double d = 0.0;
+            var sync = new object();
-
+            var d = 0.0;
-            for (int i = 0; i < y.Length; i++)
+            Parallel.For(
-            {
+                0,
-                d += y[i] * x[i];
+                y.Length,
-            }
+                () => 0.0,
+                (i, localData) => localData += y[i] * x[i],
+                localResult =>
+                {
+                    lock (sync)
+                    {
+                        d += localResult;
+                    }
+                }
+                );
 
             return d;
         }
@@ -731,23 +740,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             {
                 var d = 0.0;
                 int index;
-                Parallel.For(
+                for (var k = 0; k < j; k++)
-                    0,
+                {
-                    j,
+                    var s = 0.0;
-                    k =>
+                    int i;
-                    // for (var k = 0; k < j; k++)
+                    for (i = 0; i < k; i++)
                     {
-                        var s = 0.0;
+                        s += a[i * order + k] * a[i * order + j];
-                        int i;
+                    }
-                        for (i = 0; i < k; i++)
+                    var tmp = k * order;
-                        {
+                    index = tmp + j;
-                            s += a[i * order + k] * a[i * order + j];
+                    a[index] = s = (a[index] - s) / a[tmp + k];
-                        }
+                    d += s * s;
-                        var tmp = k * order;
+                }
-                        index = tmp + j;
-                        a[index] = s = (a[index] - s) / a[tmp + k];
-                        d += s * s;
-                    });
 
                 index = j * order + j;
                 d = a[index] - d;
@@ -756,7 +761,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                a[index] = System.Math.Sqrt(d);
+                a[index] = Math.Sqrt(d);
                 for (var k = j + 1; k < order; k++)
                 {
                     a[k * order + j] = 0.0;