Parallelized Cholesky Solve.

Bug Fix in Cholesky factorization.
Unit tests for Cholesky Solve on random matrices.

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs

@@ -742,6 +742,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
         public void CholeskyFactor(double[] a, int order)
         {
+            if (a == null)
+            {
+                throw new ArgumentNullException("a");
+            }
+
             for (var j = 0; j < order; j++)
             {
                 var d = 0.0;
@@ -793,15 +798,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <summary>
         /// Solves A*X=B for X using a previously factored A matrix.
         /// </summary>
-        /// <param name="a">The square, positive definite matrix A.</param>
+        /// <param name="a">The square, positive definite matrix A. Has to be different than <paramref name="B"/>.</param>
         /// <param name="aOrder">The number of rows and columns in A.</param>
-        /// <param name="b">The B matrix.</param>
+        /// <param name="b">The B matrix. Has to be different than <paramref name="A"/>.</param>
         /// <param name="bRows">The number of rows in the B matrix.</param>
         /// <param name="bColumns">The number of columns in the B matrix.</param>
         /// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
         public void CholeskySolveFactored(double[] a, int aOrder, double[] b, int bRows, int bColumns)
         {
-            for (int c = 0; c < bColumns; c++)
+            if (a == null)
+            {
+                throw new ArgumentNullException("a");
+            }
+
+            if (b == null)
+            {
+                throw new ArgumentNullException("b");
+            }
+
+            if (aOrder != bRows)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            if (Object.ReferenceEquals(a, b))
+            {
+                throw new ArgumentException(Resources.ArgumentReferenceDifferent);
+            }
+
+            CommonParallel.For(0, bColumns, c =>
             {
                 int cindex = c * aOrder;
 
@@ -809,7 +834,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 double sum;
                 for (int i = 0; i < aOrder; i++)
                 {
-                    sum = b[c * aOrder + i];
+                    sum = b[cindex + i];
                     for (int k = i - 1; k >= 0; k--)
                     {
                         sum -= a[k * aOrder + i] * b[cindex + k];
@@ -828,7 +853,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                     }
                     b[cindex + i] = sum / a[iindex + i];
                 }
-            }
+            });
         }
 
         public void QRFactor(double[] r, double[] q)

src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs

@@ -128,7 +128,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment.");
             }
 
-            var dresult = input as DenseMatrix;
+            var dresult = result as DenseMatrix;
             if (dresult == null)
             {
                 throw new NotImplementedException("Can only do Cholesky factorization for dense matrices at the moment.");
@@ -195,7 +195,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment.");
             }
 
-            var dresult = input as DenseVector;
+            var dresult = result as DenseVector;
             if (dresult == null)
             {
                 throw new NotImplementedException("Can only do Cholesky factorization for dense vectors at the moment.");

src/Numerics/Properties/Resources.Designer.cs

@@ -384,6 +384,15 @@ namespace MathNet.Numerics.Properties {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized string similar to Arguments must be different objects..
+        /// </summary>
+        internal static string ArgumentReferenceDifferent {
+            get {
+                return ResourceManager.GetString("ArgumentReferenceDifferent", resourceCulture);
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized string similar to Array must have exactly one dimension (and not be null)..
         /// </summary>

src/Numerics/Properties/Resources.resx

@@ -112,10 +112,10 @@
     <value>2.0</value>
   </resheader>
   <resheader name="reader">
-    <value>System.Resources.ResXResourceReader, System.Windows.Forms, Version=2.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089</value>
+    <value>System.Resources.ResXResourceReader, System.Windows.Forms, Version=4.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089</value>
   </resheader>
   <resheader name="writer">
-    <value>System.Resources.ResXResourceWriter, System.Windows.Forms, Version=2.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089</value>
+    <value>System.Resources.ResXResourceWriter, System.Windows.Forms, Version=4.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089</value>
   </resheader>
   <data name="ArgumentHistogramContainsNot" xml:space="preserve">
     <value>The histogram does not contains the value.</value>
@@ -297,4 +297,7 @@
   <data name="UndefinedMoment" xml:space="preserve">
     <value>The moment of the distribution is undefined.</value>
   </data>
+  <data name="ArgumentReferenceDifferent" xml:space="preserve">
+    <value>Arguments must be different objects.</value>
+  </data>
 </root>

src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs

@@ -103,20 +103,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanFactorizeRandomMatrix(int order)
         {
-            // Fill a matrix with standard random numbers.
-            var normal = new Distributions.Normal();
-            normal.RandomSource = new Random.MersenneTwister(1);
-            var A = new DenseMatrix(order);
-            for (int i = 0; i < order; i++)
-            {
-                for (int j = 0; j < order; j++)
-                {
-                    A[i, j] = normal.Sample();
-                }
-            }
-
-            // Generate a matrix which is positive definite.
-            var X = A.Transpose() * A;
+            var X = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order);
             var chol = X.Cholesky();
             var C = chol.Factor;
 
@@ -139,7 +126,178 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
             {
                 for (int j = 0; j < XfromC.ColumnCount; j++)
                 {
-                    Assert.AreApproximatelyEqual(X[i,j], XfromC[i, j], 1.0e-13);
+                    Assert.AreApproximatelyEqual(X[i,j], XfromC[i, j], 1.0e-11);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomVector(int order)
+        {
+            var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order);
+            var ACopy = A.Clone();
+            var chol = A.Cholesky();
+            var b = MatrixLoader.GenerateRandomVector(order);
+            var x = chol.Solve(b);
+
+            Assert.AreEqual(b.Count, x.Count);
+
+            var bReconstruct = A * x;
+
+            // Check the reconstruction.
+            for (int i = 0; i < order; i++)
+            {
+                Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (int i = 0; i < A.RowCount; i++)
+            {
+                for (int j = 0; j < A.ColumnCount; j++)
+                {
+                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1,1)]
+        [Row(2,4)]
+        [Row(5,8)]
+        [Row(10,3)]
+        [Row(50,10)]
+        [Row(100,100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrix(int row, int col)
+        {
+            var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(row);
+            var ACopy = A.Clone();
+            var chol = A.Cholesky();
+            var B = MatrixLoader.GenerateRandomMatrix(row, col);
+            var X = chol.Solve(B);
+
+            Assert.AreEqual(B.RowCount, X.RowCount);
+            Assert.AreEqual(B.ColumnCount, X.ColumnCount);
+
+            var BReconstruct = A * X;
+
+            // Check the reconstruction.
+            for (int i = 0; i < B.RowCount; i++)
+            {
+                for (int j = 0; j < B.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (int i = 0; i < A.RowCount; i++)
+            {
+                for (int j = 0; j < A.ColumnCount; j++)
+                {
+                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomVectorWhenResultVectorGiven(int order)
+        {
+            var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order);
+            var ACopy = A.Clone();
+            var chol = A.Cholesky();
+            var b = MatrixLoader.GenerateRandomVector(order);
+            var bCopy = b.Clone();
+            var x = new DenseVector(order);
+            chol.Solve(b, x);
+
+            Assert.AreEqual(b.Count, x.Count);
+
+            var bReconstruct = A * x;
+
+            // Check the reconstruction.
+            for (int i = 0; i < order; i++)
+            {
+                Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (int i = 0; i < A.RowCount; i++)
+            {
+                for (int j = 0; j < A.ColumnCount; j++)
+                {
+                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                }
+            }
+
+            // Make sure b didn't change.
+            for (int i = 0; i < order; i++)
+            {
+                Assert.AreEqual(bCopy[i], b[i]);
+            }
+        }
+
+        [Test]
+        [Row(1, 1)]
+        [Row(2, 4)]
+        [Row(5, 8)]
+        [Row(10, 3)]
+        [Row(50, 10)]
+        [Row(100, 100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int col)
+        {
+            var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(row);
+            var ACopy = A.Clone();
+            var chol = A.Cholesky();
+            var B = MatrixLoader.GenerateRandomMatrix(row, col);
+            var BCopy = B.Clone();
+            var X = new DenseMatrix(row, col);
+            chol.Solve(B, X);
+
+            Assert.AreEqual(B.RowCount, X.RowCount);
+            Assert.AreEqual(B.ColumnCount, X.ColumnCount);
+
+            var BReconstruct = A * X;
+
+            // Check the reconstruction.
+            for (int i = 0; i < B.RowCount; i++)
+            {
+                for (int j = 0; j < B.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (int i = 0; i < A.RowCount; i++)
+            {
+                for (int j = 0; j < A.ColumnCount; j++)
+                {
+                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                }
+            }
+
+            // Make sure B didn't change.
+            for (int i = 0; i < B.RowCount; i++)
+            {
+                for (int j = 0; j < B.ColumnCount; j++)
+                {
+                    Assert.AreEqual(BCopy[i, j], B[i, j]);
                 }
             }
         }

src/UnitTests/LinearAlgebraTests/Double/MatrixLoader.cs

@@ -63,5 +63,55 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             }
         }
 
+        public static Matrix GenerateRandomMatrix(int row, int col)
+        {
+            // Fill a matrix with standard random numbers.
+            var normal = new Distributions.Normal();
+            normal.RandomSource = new Random.MersenneTwister(1);
+            var A = new DenseMatrix(row, col);
+            for (int i = 0; i < row; i++)
+            {
+                for (int j = 0; j < col; j++)
+                {
+                    A[i, j] = normal.Sample();
+                }
+            }
+
+            // Generate a matrix which is positive definite.
+            return A;
+        }
+
+        public static Matrix GenerateRandomPositiveDefiniteMatrix(int order)
+        {
+            // Fill a matrix with standard random numbers.
+            var normal = new Distributions.Normal();
+            normal.RandomSource = new Random.MersenneTwister(1);
+            var A = new DenseMatrix(order);
+            for (int i = 0; i < order; i++)
+            {
+                for (int j = 0; j < order; j++)
+                {
+                    A[i, j] = normal.Sample();
+                }
+            }
+
+            // Generate a matrix which is positive definite.
+            return A.Transpose() * A;
+        }
+
+        public static Vector GenerateRandomVector(int order)
+        {
+            // Fill a matrix with standard random numbers.
+            var normal = new Distributions.Normal();
+            normal.RandomSource = new Random.MersenneTwister(1);
+            var v = new DenseVector(order);
+            for (int i = 0; i < order; i++)
+            {
+                v[i] = normal.Sample();
+            }
+
+            // Generate a matrix which is positive definite.
+            return v;
+        }
     }
 }