linear algebra: modified cholesky to run inplace (removed copying to another array)

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs

@@ -727,39 +727,41 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
         public void CholeskyFactor(double[] a, int order)
         {
-            double[] factor = new double[a.Length];
+            for (var j = 0; j < order; j++)
-
-            for (int j = 0; j < order; j++)
             {
-                double d = 0.0;
+                var d = 0.0;
                 int index;
-                for (int k = 0; k < j; k++)
+                Parallel.For(
-                {
+                    0,
-                    double s = 0.0;
+                    j,
-                    int i;
+                    k =>
-                    for (i = 0; i < k; i++)
+                    // for (var k = 0; k < j; k++)
                     {
-                        s += factor[i * order + k] * factor[i * order + j];
+                        var s = 0.0;
-                    }
+                        int i;
-                    int tmp = k * order;
+                        for (i = 0; i < k; i++)
-                    index = tmp + j;
+                        {
-                    factor[index] = s = (a[index] - s) / factor[tmp + k];
+                            s += a[i * order + k] * a[i * order + j];
-                    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;
                 if (d <= 0.0)
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
-                factor[index] = System.Math.Sqrt(d);
+
-                for (int k = j + 1; k < order; k++)
+                a[index] = System.Math.Sqrt(d);
+                for (var k = j + 1; k < order; k++)
                 {
-                    factor[k * order + j] = 0.0;
+                    a[k * order + j] = 0.0;
                 }
             }
-
-            Buffer.BlockCopy(factor, 0, a, 0, factor.Length * Constants.SizeOfDouble);
         }
 
         public void CholeskySolve(int columnsOfB, double[] a, double[] b)

src/UnitTests/LinearAlgebraTests/Double/LinearAlgebraProviderTests.cs

@@ -2,6 +2,7 @@
 {
     using System;
     using Algorithms.LinearAlgebra;
+    using LinearAlgebra.Double;
     using MbUnit.Framework;
 
     [TestFixture]
@@ -176,10 +177,27 @@
 
         }
 
-        [Test, Ignore]
+        [Test]
-        public void CanComputeCholeskyFactor(double[] a, int order)
+        public void CanComputeCholeskyFactor()
         {
-
+            var matrix = new double[] { 1, 1, 1, 1, 1, 5, 5, 5, 1, 5, 14, 14, 1, 5, 14, 15 };
+            Provider.CholeskyFactor(matrix, 4);
+            Assert.AreEqual(matrix[0], 1);
+            Assert.AreEqual(matrix[1], 1);
+            Assert.AreEqual(matrix[2], 1);
+            Assert.AreEqual(matrix[3], 1);
+            Assert.AreEqual(matrix[4], 0);
+            Assert.AreEqual(matrix[5], 2);
+            Assert.AreEqual(matrix[6], 2);
+            Assert.AreEqual(matrix[7], 2);
+            Assert.AreEqual(matrix[8], 0);
+            Assert.AreEqual(matrix[9], 0);
+            Assert.AreEqual(matrix[10], 3);
+            Assert.AreEqual(matrix[11], 3);
+            Assert.AreEqual(matrix[12], 0);
+            Assert.AreEqual(matrix[13], 0);
+            Assert.AreEqual(matrix[14], 0);
+            Assert.AreEqual(matrix[15], 1);
         }
 
         [Test, Ignore]