factor: Merged Andriy's factorization code for user defined matrices

src/MathNet.Numerics.5.0.ReSharper

@@ -863,6 +863,7 @@ Cholesky</UserWords>
       <Abbreviation Text="FFT" />
       <Abbreviation Text="LU" />
       <Abbreviation Text="DC" />
+      <Abbreviation Text="VT" />
     </Naming2>
   </CodeStyleSettings>
   <SharedSolutionTemplateManager>

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs

@@ -1361,7 +1361,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <summary>
         /// Perform calculation of Q or R
         /// </summary>
-        /// <param name="work">Work arrat</param>
+        /// <param name="work">Work array</param>
         /// <param name="workIndex">Index of colunn in work array</param>
         /// <param name="a">Q or R matrices</param>
         /// <param name="rowCount">The number of rows</param>
@@ -3053,16 +3053,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         }
 
         /// <summary>
-        /// Сonstruct givens plane rotation
+        /// Given the Cartesian coordinates (da, db) of a point p, these fucntion return the parameters da, db, c, and s 
+        /// associated with the Givens rotation that zeros the y-coordinate of the point.
         /// </summary>
-        /// <param name="da"></param>
-        /// <param name="db"></param>
-        /// <param name="c"></param>
-        /// <param name="s"></param>
+        /// <param name="da">Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation</param>
+        /// <param name="db">Provides the y-coordinate of the point p. On exit contains the parameter z associated with the Givens rotation</param>
+        /// <param name="c">Contains the parameter c associated with the Givens rotation</param>
+        /// <param name="s">Contains the parameter s associated with the Givens rotation</param>
+        /// <remarks>This is equivalent to the DROTG LAPACK routine.</remarks>
         private static void Drotg(ref double da, ref double db, ref double c, ref double s)
         {
-            // Сonstruct givens plane rotation.
-            // jack dongarra, linpack, 3/11/78.
             double r, z;
 
             var roe = db;
@@ -3107,7 +3107,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             da = r;
             db = z;
-            return;
         }
 
         /// <summary>

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

@@ -56,16 +56,27 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 return new DenseCholesky(dense);
             }
 
-            throw new NotImplementedException();
+            return new UserCholesky(matrix);
         }
 
         /// <summary>
-        /// Gets or sets the lower triangular form of the Cholesky matrix.
+        /// Gets or sets the lower triangular form of the Cholesky matrix
         /// </summary>
-        public virtual Matrix Factor
+        protected Matrix CholeskyFactor
         {
             get;
-            protected set;
+            set;
+        }
+
+        /// <summary>
+        /// Gets the lower triangular form of the Cholesky matrix.
+        /// </summary>
+        public virtual Matrix Factor
+        {
+            get
+            {
+                return CholeskyFactor.Clone();
+            }
         }
 
         /// <summary>
@@ -76,9 +87,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             get
             {
                 var det = 1.0;
-                for (var j = 0; j < Factor.RowCount; j++)
+                for (var j = 0; j < CholeskyFactor.RowCount; j++)
                 {
-                    det *= Factor[j, j] * Factor[j, j];
+                    det *= CholeskyFactor[j, j] * CholeskyFactor[j, j];
                 }
 
                 return det;
@@ -93,9 +104,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             get
             {
                 var det = 0.0;
-                for (var j = 0; j < Factor.RowCount; j++)
+                for (var j = 0; j < CholeskyFactor.RowCount; j++)
                 {
-                    det += 2.0 * Math.Log(Factor[j, j]);
+                    det += 2.0 * Math.Log(CholeskyFactor[j, j]);
                 }
 
                 return det;

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

@@ -67,7 +67,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
             var factor = (DenseMatrix)matrix.Clone();
             Control.LinearAlgebraProvider.CholeskyFactor(factor.Data, factor.RowCount);
-            Factor = factor;
+            CholeskyFactor = factor;
         }
 
         /// <summary>
@@ -99,7 +99,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
             }
 
-            if (input.RowCount != Factor.RowCount)
+            if (input.RowCount != CholeskyFactor.RowCount)
             {
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
@@ -120,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble);
 
             // Cholesky solve by overwriting result.
-            var dfactor = (DenseMatrix)Factor;
+            var dfactor = (DenseMatrix)CholeskyFactor;
             Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount);
         }
 
@@ -148,7 +148,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-            if (input.Count != Factor.RowCount)
+            if (input.Count != CholeskyFactor.RowCount)
             {
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
@@ -169,7 +169,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble);
 
             // Cholesky solve by overwriting result.
-            var dfactor = (DenseMatrix)Factor;
+            var dfactor = (DenseMatrix)CholeskyFactor;
             Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1);
         }
     }

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

@@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 return new DenseLU(dense);
             }
 
-            throw new NotImplementedException();
+            return new UserLU(matrix);
         }
 
         /// <summary>

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

@@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 return new DenseQR(dense);
             }
 
-            throw new NotImplementedException();
+            return new UserQR(matrix);
         }
 
         /// <summary>

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

@@ -123,7 +123,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 return new DenseSvd(dense, computeVectors);
             }
 
-            throw new NotImplementedException();
+            return new UserSvd(matrix, computeVectors);
         }
 
         /// <summary>

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

@@ -0,0 +1,222 @@
+// <copyright file="UserCholesky.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebra.Double.Factorization
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
+    /// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
+    /// is an lower triangular matrix L so that A = L*L'.</para>
+    /// </summary>
+    /// <remarks>
+    /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
+    /// or positive definite, the constructor will throw an exception.
+    /// </remarks>
+    public class UserCholesky : Cholesky
+    {
+        /// <summary>
+        /// Initializes a new instance of the <see cref="UserCholesky"/> class. This object will compute the
+        /// Cholesky factorization when the constructor is called and cache it's factorization.
+        /// </summary>
+        /// <param name="matrix">The matrix to factor.</param>
+        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
+        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
+        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not positive definite.</exception>
+        public UserCholesky(Matrix matrix)
+        {
+            if (matrix == null)
+            {
+                throw new ArgumentNullException("matrix");
+            }
+
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare);
+            }
+
+            // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
+            CholeskyFactor = matrix.Clone();
+            for (var j = 0; j < CholeskyFactor.RowCount; j++)
+            {
+                var d = 0.0;
+                for (var k = 0; k < j; k++)
+                {
+                    var s = 0.0;
+                    for (var i = 0; i < k; i++)
+                    {
+                        s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i);
+                    }
+
+                    s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k);
+                    CholeskyFactor.At(j, k, s);
+                    d += s * s;
+                }
+
+                d = matrix.At(j, j) - d;
+                if (d <= 0.0)
+                {
+                    throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
+                }
+
+                CholeskyFactor.At(j, j, Math.Sqrt(d));
+                for (var k = j + 1; k < CholeskyFactor.RowCount; k++)
+                {
+                    CholeskyFactor.At(j, k, 0.0);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Solves a system of linear equations, <b>AX = B</b>, with A Cholesky factorized.
+        /// </summary>
+        /// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
+        public override void Solve(Matrix input, Matrix result)
+        {
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            // Check for proper dimensions.
+            if (result.RowCount != input.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+            }
+
+            if (result.ColumnCount != input.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+            }
+
+            if (input.RowCount != CholeskyFactor.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            input.CopyTo(result);
+            var order = CholeskyFactor.RowCount;
+
+            for (var c = 0; c < result.ColumnCount; c++)
+            {
+                // Solve L*Y = B;
+                double sum;
+                for (var i = 0; i < order; i++)
+                {
+                    sum = result.At(i, c);
+                    for (var k = i - 1; k >= 0; k--)
+                    {
+                        sum -= CholeskyFactor.At(i, k) * result.At(k, c);
+                    }
+
+                    result.At(i, c, sum / CholeskyFactor.At(i, i));
+                }
+
+                // Solve L'*X = Y;
+                for (var i = order - 1; i >= 0; i--)
+                {
+                    sum = result.At(i, c);
+                    for (var k = i + 1; k < order; k++)
+                    {
+                        sum -= CholeskyFactor.At(k, i) * result.At(k, c);
+                    }
+
+                    result.At(i, c, sum / CholeskyFactor.At(i, i));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Solves a system of linear equations, <b>Ax = b</b>, with A Cholesky factorized.
+        /// </summary>
+        /// <param name="input">The right hand side vector, <b>b</b>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>x</b>.</param>
+        public override void Solve(Vector input, Vector result)
+        {
+            // Check for proper arguments.
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            // Check for proper dimensions.
+            if (input.Count != result.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (input.Count != CholeskyFactor.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            input.CopyTo(result);
+            var order = CholeskyFactor.RowCount;
+
+            // Solve L*Y = B;
+            double sum;
+            for (var i = 0; i < order; i++)
+            {
+                sum = result[i];
+                for (var k = i - 1; k >= 0; k--)
+                {
+                    sum -= CholeskyFactor.At(i, k) * result[k];
+                }
+
+                result[i] = sum / CholeskyFactor.At(i, i);
+            }
+
+            // Solve L'*X = Y;
+            for (var i = order - 1; i >= 0; i--)
+            {
+                sum = result[i];
+                for (var k = i + 1; k < order; k++)
+                {
+                    sum -= CholeskyFactor.At(k, i) * result[k];
+                }
+
+                result[i] = sum / CholeskyFactor.At(i, i);
+            }
+        }
+    }
+}

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

@@ -0,0 +1,300 @@
+// <copyright file="UserLU.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebra.Double.Factorization
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// <para>A class which encapsulates the functionality of an LU factorization.</para>
+    /// <para>For a matrix A, the LU factorization is a pair of lower triangular matrix L and
+    /// upper triangular matrix U so that A = L*U.</para>
+    /// </summary>
+    /// <remarks>
+    /// The computation of the LU factorization is done at construction time.
+    /// </remarks>
+    public class UserLU : LU
+    {
+        /// <summary>
+        /// Initializes a new instance of the <see cref="UserLU"/> class. This object will compute the
+        /// LU factorization when the constructor is called and cache it's factorization.
+        /// </summary>
+        /// <param name="matrix">The matrix to factor.</param>
+        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
+        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
+        public UserLU(Matrix matrix)
+        {
+            if (matrix == null)
+            {
+                throw new ArgumentNullException("matrix");
+            }
+
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare);
+            }
+
+            // Create an array for the pivot indices.
+            var order = matrix.RowCount;
+            Factors = matrix.Clone();
+            Pivots = new int[order];
+            
+            // Initialize the pivot matrix to the identity permutation.
+            for (var i = 0; i < order; i++)
+            {
+                Pivots[i] = i;
+            }
+
+            var vectorLUcolj = new double[order];
+            for (var j = 0; j < order; j++)
+            {
+                // Make a copy of the j-th column to localize references.
+                for (var i = 0; i < order; i++)
+                {
+                    vectorLUcolj[i] = Factors.At(i, j);
+                }
+
+                // Apply previous transformations.
+                for (var i = 0; i < order; i++)
+                {
+                    var kmax = Math.Min(i, j);
+                    var s = 0.0;
+                    for (var k = 0; k < kmax; k++)
+                    {
+                        s += Factors.At(i, k) * vectorLUcolj[k];
+                    }
+
+                    vectorLUcolj[i] -= s;
+                    Factors.At(i, j, vectorLUcolj[i]);
+                }
+
+                // Find pivot and exchange if necessary.
+                var p = j;
+                for (var i = j + 1; i < order; i++)
+                {
+                    if (Math.Abs(vectorLUcolj[i]) > Math.Abs(vectorLUcolj[p]))
+                    {
+                        p = i;
+                    }
+                }
+
+                if (p != j)
+                {
+                    for (var k = 0; k < order; k++)
+                    {
+                        var temp = Factors.At(p, k);
+                        Factors.At(p, k, Factors.At(j, k));
+                        Factors.At(j, k, temp);
+                    }
+
+                    Pivots[j] = p;
+                }
+
+                // Compute multipliers.
+                if (j < order & Factors.At(j, j) != 0.0)
+                {
+                    for (var i = j + 1; i < order; i++)
+                    {
+                        Factors.At(i, j, (Factors.At(i, j) / Factors.At(j, j)));
+                    }
+                }
+            }
+        }
+
+        /// <summary>
+        /// Solves a system of linear equations, <c>AX = B</c>, with A LU factorized.
+        /// </summary>
+        /// <param name="input">The right hand side <see cref="Matrix"/>, <c>B</c>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <c>X</c>.</param>
+        public override void Solve(Matrix input, Matrix result)
+        {
+            // Check for proper arguments.
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            // Check for proper dimensions.
+            if (result.RowCount != input.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+            }
+
+            if (result.ColumnCount != input.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+            }
+
+            if (input.RowCount != Factors.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            // Copy the contents of input to result.
+            input.CopyTo(result);
+            for (var i = 0; i < Pivots.Length; i++)
+            {
+                if (Pivots[i] == i)
+                {
+                    continue;
+                }
+
+                var p = Pivots[i];
+                for (var j = 0; j < result.ColumnCount; j++)
+                {
+                    var temp = result.At(p, j);
+                    result.At(p, j, result.At(i, j));
+                    result.At(i, j, temp);
+                }
+            }
+
+            var order = Factors.RowCount;
+            
+            // Solve L*Y = P*B
+            for (var k = 0; k < order; k++)
+            {
+                for (var i = k + 1; i < order; i++)
+                {
+                    for (var j = 0; j < result.ColumnCount; j++)
+                    {
+                        var temp = result.At(k, j) * Factors.At(i, k);
+                        result.At(i, j, result.At(i, j) - temp);
+                    }
+                }
+            }
+
+            // Solve U*X = Y;
+            for (var k = order - 1; k >= 0; k--)
+            {
+                for (var j = 0; j < result.ColumnCount; j++)
+                {
+                    result.At(k, j, (result.At(k, j) / Factors.At(k, k)));
+                }
+
+                for (var i = 0; i < k; i++)
+                {
+                    for (var j = 0; j < result.ColumnCount; j++)
+                    {
+                        var temp = result.At(k, j) * Factors.At(i, k);
+                        result.At(i, j, result.At(i, j) - temp);
+                    }
+                }
+            }
+        }
+
+        /// <summary>
+        /// Solves a system of linear equations, <c>Ax = b</c>, with A LU factorized.
+        /// </summary>
+        /// <param name="input">The right hand side vector, <c>b</c>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <c>x</c>.</param>
+        public override void Solve(Vector input, Vector result)
+        {
+            // Check for proper arguments.
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            // Check for proper dimensions.
+            if (input.Count != result.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (input.Count != Factors.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            // Copy the contents of input to result.
+            input.CopyTo(result);
+            for (var i = 0; i < Pivots.Length; i++)
+            {
+                if (Pivots[i] == i)
+                {
+                    continue;
+                }
+
+                var p = Pivots[i];
+                var temp = result[p];
+                result[p] = result[i];
+                result[i] = temp;
+            }
+            
+            var order = Factors.RowCount;
+
+            // Solve L*Y = P*B
+            for (var k = 0; k < order; k++)
+            {
+                for (var i = k + 1; i < order; i++)
+                {
+                    result[i] -= result[k] * Factors.At(i, k);
+                }
+            }
+
+            // Solve U*X = Y;
+            for (var k = order - 1; k >= 0; k--)
+            {
+                result[k] /= Factors.At(k, k);
+                for (var i = 0; i < k; i++)
+                {
+                    result[i] -= result[k] * Factors.At(i, k);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition.
+        /// </summary>
+        /// <returns>The inverse of this matrix.</returns>
+        public override Matrix Inverse()
+        {
+            var order = Factors.RowCount;
+            var inverse = Factors.CreateMatrix(order, order);
+            for (var i = 0; i < order; i++)
+            {
+                inverse.At(i, i, 1.0);
+            }
+
+            return Solve(inverse);
+        }
+    }
+}

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

@@ -0,0 +1,332 @@
+// <copyright file="UserQR.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebra.Double.Factorization
+{
+    using System;
+    using System.Linq;
+    using Properties;
+
+    /// <summary>
+    /// <para>A class which encapsulates the functionality of the QR decomposition.</para>
+    /// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal matrix 
+    /// (its columns are orthogonal unit vectors meaning QTQ = I) and R is an upper triangular matrix 
+    /// (also called right triangular matrix).</para>
+    /// </summary>
+    /// <remarks>
+    /// The computation of the QR decomposition is done at construction time by Householder transformation.
+    /// </remarks>
+    public class UserQR : QR
+    {
+        /// <summary>
+        /// Initializes a new instance of the <see cref="UserQR"/> class. This object will compute the
+        /// QR factorization when the constructor is called and cache it's factorization.
+        /// </summary>
+        /// <param name="matrix">The matrix to factor.</param>
+        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
+        public UserQR(Matrix matrix)
+        {
+            if (matrix == null)
+            {
+                throw new ArgumentNullException("matrix");
+            }
+
+            if (matrix.RowCount < matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            MatrixR = matrix.Clone();
+            MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount);
+
+            for (var i = 0; i < matrix.RowCount; i++)
+            {
+                MatrixQ.At(i, i, 1.0);
+            }
+
+            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
+            var u = new double[minmn][];
+            for (var i = 0; i < minmn; i++)
+            {
+                u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
+                ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
+            }
+
+            for (var i = minmn - 1; i >= 0; i--)
+            {
+                ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
+            }
+        }
+
+        /// <summary>
+        /// Generate column from initial matrix to work array
+        /// </summary>
+        /// <param name="a">Initial matrix</param>
+        /// <param name="rowStart">The firts row</param>
+        /// <param name="rowEnd">The last row</param>
+        /// <param name="column">Column index</param>
+        /// <returns>Generated vector</returns>
+        private static double[] GenerateColumn(Matrix a, int rowStart, int rowEnd, int column)
+        {
+            var ru = rowEnd - rowStart + 1;
+            var u = new double[ru];
+
+            for (var i = rowStart; i <= rowEnd; i++)
+            {
+                u[i - rowStart] = a.At(i, rowStart);
+                a.At(i, rowStart, 0.0);
+            }
+
+            var norm = u.Sum(t => t * t);
+            norm = Math.Sqrt(norm);
+
+            if (rowStart == rowEnd || norm == 0)
+            {
+                a.At(rowStart, column, -u[0]);
+                u[0] = Math.Sqrt(2.0);
+                return u;
+            }
+
+            var scale = 1.0 / norm;
+            if (u[0] < 0.0)
+            {
+                scale *= -1.0;
+            }
+
+            a.At(rowStart, column, -1.0 / scale);
+
+            for (var i = 0; i < ru; i++)
+            {
+                u[i] *= scale;
+            }
+
+            u[0] += 1.0;
+            var s = Math.Sqrt(1.0 / u[0]);
+
+            for (var i = 0; i < ru; i++)
+            {
+                u[i] *= s;
+            }
+
+            return u;
+        }
+
+        /// <summary>
+        /// Perform calculation of Q or R
+        /// </summary>
+        /// <param name="u">Work array</param>
+        /// <param name="a">Q or R matrices</param>
+        /// <param name="rowStart">The first row</param>
+        /// <param name="rowEnd">The last row</param>
+        /// <param name="columnStart">The first column</param>
+        /// <param name="columnEnd">The last column</param>
+        private static void ComputeQR(double[] u, Matrix a, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        {
+            if (rowEnd < rowStart || columnEnd < columnStart)
+            {
+                return;
+            }
+
+            var v = new double[columnEnd - columnStart + 1];
+            for (var j = columnStart; j <= columnEnd; j++)
+            {
+                v[j - columnStart] = 0.0;
+            }
+
+            for (var i = rowStart; i <= rowEnd; i++)
+            {
+                for (var j = columnStart; j <= columnEnd; j++)
+                {
+                    v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
+                }
+            }
+
+            for (var i = rowStart; i <= rowEnd; i++)
+            {
+                for (var j = columnStart; j <= columnEnd; j++)
+                {
+                    a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart]));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Solves a system of linear equations, <b>AX = B</b>, with A QR factorized.
+        /// </summary>
+        /// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
+        public override void Solve(Matrix input, Matrix result)
+        {
+            // Check for proper arguments.
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            // The solution X should have the same number of columns as B
+            if (input.ColumnCount != result.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+            }
+
+            // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
+            if (MatrixR.RowCount != input.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+            }
+
+            // The solution X row dimension is equal to the column dimension of A
+            if (MatrixR.ColumnCount != result.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+            }
+
+            var inputCopy = input.Clone();
+
+            // Compute Y = transpose(Q)*B
+            var bn = inputCopy.ColumnCount;
+            var column = new double[MatrixR.RowCount];
+            for (var j = 0; j < bn; j++)
+            {
+                for (var k = 0; k < MatrixR.RowCount; k++)
+                {
+                    column[k] = inputCopy.At(k, j);
+                }
+
+                for (var i = 0; i < MatrixR.RowCount; i++)
+                {
+                    double s = 0;
+                    for (var k = 0; k < MatrixR.RowCount; k++)
+                    {
+                        s += MatrixQ.At(k, i) * column[k];
+                    }
+
+                    inputCopy.At(i, j, s);
+                }
+            }
+
+            // Solve R*X = Y;
+            for (var k = MatrixR.ColumnCount - 1; k >= 0; k--)
+            {
+                for (var j = 0; j < bn; j++)
+                {
+                    inputCopy.At(k, j, inputCopy.At(k, j) / MatrixR.At(k, k));
+                }
+
+                for (var i = 0; i < k; i++)
+                {
+                    for (var j = 0; j < bn; j++)
+                    {
+                        inputCopy.At(i, j, inputCopy.At(i, j) - (inputCopy.At(k, j) * MatrixR.At(i, k)));
+                    }
+                }
+            }
+
+            for (var i = 0; i < MatrixR.ColumnCount; i++)
+            {
+                for (var j = 0; j < inputCopy.ColumnCount; j++)
+                {
+                    result.At(i, j, inputCopy.At(i, j));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Solves a system of linear equations, <b>Ax = b</b>, with A QR factorized.
+        /// </summary>
+        /// <param name="input">The right hand side vector, <b>b</b>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>x</b>.</param>
+        public override void Solve(Vector input, Vector result)
+        {
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            // Ax=b where A is an m x n matrix
+            // Check that b is a column vector with m entries
+            if (MatrixR.RowCount != input.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            // Check that x is a column vector with n entries
+            if (MatrixR.ColumnCount != result.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            var inputCopy = input.Clone();
+
+            // Compute Y = transpose(Q)*B
+            var column = new double[MatrixR.RowCount];
+            for (var k = 0; k < MatrixR.RowCount; k++)
+            {
+                column[k] = inputCopy[k];
+            }
+
+            for (var i = 0; i < MatrixR.RowCount; i++)
+            {
+                double s = 0;
+                for (var k = 0; k < MatrixR.RowCount; k++)
+                {
+                    s += MatrixQ.At(k, i) * column[k];
+                }
+
+                inputCopy[i] = s;
+            }
+
+            // Solve R*X = Y;
+            for (var k = MatrixR.ColumnCount - 1; k >= 0; k--)
+            {
+                inputCopy[k] /= MatrixR.At(k, k);
+                for (var i = 0; i < k; i++)
+                {
+                    inputCopy[i] -= inputCopy[k] * MatrixR.At(i, k);
+                }
+            }
+
+            for (var i = 0; i < MatrixR.ColumnCount; i++)
+            {
+                result[i] = inputCopy[i];
+            }
+        }
+    }
+}

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

@@ -0,0 +1,923 @@
+// <copyright file="UserSvd.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebra.Double.Factorization
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// <para>A class which encapsulates the functionality of the singular value decomposition (SVD) for <see cref="Matrix"/>.</para>
+    /// <para>Suppose M is an m-by-n matrix whose entries are real numbers. 
+    /// Then there exists a factorization of the form M = UΣVT where:
+    /// - U is an m-by-m unitary matrix;
+    /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
+    /// - VT denotes transpose of V, an n-by-n unitary matrix; 
+    /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal 
+    /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined 
+    /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
+    /// </summary>
+    /// <remarks>
+    /// The computation of the singular value decomposition is done at construction time.
+    /// </remarks>
+    public class UserSvd : Svd
+    {
+        /// <summary>
+        /// Initializes a new instance of the <see cref="UserSvd"/> class. This object will compute the
+        /// the singular value decomposition when the constructor is called and cache it's decomposition.
+        /// </summary>
+        /// <param name="matrix">The matrix to factor.</param>
+        /// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
+        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
+        /// <exception cref="ArgumentException">If SVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
+        public UserSvd(Matrix matrix, bool computeVectors)
+        {
+            if (matrix == null)
+            {
+                throw new ArgumentNullException("matrix");
+            }
+
+            ComputeVectors = computeVectors;
+            var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount);
+            var matrixCopy = matrix.Clone();
+
+            VectorS = matrixCopy.CreateVector(nm);
+            MatrixU = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount);
+            MatrixVT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount);
+
+            const int Maxiter = 1000;
+            var e = new double[matrixCopy.ColumnCount];
+            var work = new double[matrixCopy.RowCount];
+
+            int i, j;
+            int l, lp1;
+            var cs = 0.0;
+            var sn = 0.0;
+            double t;
+
+            var ncu = matrixCopy.RowCount;
+
+            // Reduce matrixCopy to bidiagonal form, storing the diagonal elements
+            // In s and the super-diagonal elements in e.
+            var nct = Math.Min(matrixCopy.RowCount - 1, matrixCopy.ColumnCount);
+            var nrt = Math.Max(0, Math.Min(matrixCopy.ColumnCount - 2, matrixCopy.RowCount));
+            var lu = Math.Max(nct, nrt);
+            for (l = 0; l < lu; l++)
+            {
+                lp1 = l + 1;
+                if (l < nct)
+                {
+                    // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l].
+                    var xnorm = Dnrm2Column(matrixCopy, matrixCopy.RowCount, l, l);
+                    VectorS[l] = xnorm;
+                    if (VectorS[l] != 0.0)
+                    {
+                        if (matrixCopy.At(l, l) != 0.0)
+                        {
+                            VectorS[l] = Dsign(VectorS[l], matrixCopy.At(l, l));
+                        }
+
+                        DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0 / VectorS[l]);
+                        matrixCopy.At(l, l, (1.0 + matrixCopy.At(l, l)));
+                    }
+
+                    VectorS[l] = -VectorS[l];
+                }
+
+                for (j = lp1; j < matrixCopy.ColumnCount; j++)
+                {
+                    if (l < nct)
+                    {
+                        if (VectorS[l] != 0.0)
+                        {
+                            // Apply the transformation.
+                            t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l);
+                            for (var ii = l; ii < matrixCopy.RowCount; ii++)
+                            {
+                                matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l)));
+                            }
+                        }
+                    }
+
+                    // Place the l-th row of matrixCopy into  e for the
+                    // Subsequent calculation of the row transformation.
+                    e[j] = matrixCopy.At(l, j);
+                }
+
+                if (ComputeVectors && l < nct)
+                {
+                    // Place the transformation in u for subsequent back multiplication.
+                    for (i = l; i < matrixCopy.RowCount; i++)
+                    {
+                        MatrixU.At(i, l, matrixCopy.At(i, l));
+                    }
+                }
+
+                if (l >= nrt)
+                {
+                    continue;
+                }
+
+                // Compute the l-th row transformation and place the l-th super-diagonal in e(l).
+                var enorm = Dnrm2Vector(e, lp1);
+                e[l] = enorm;
+                if (e[l] != 0.0)
+                {
+                    if (e[lp1] != 0.0)
+                    {
+                        e[l] = Dsign(e[l], e[lp1]);
+                    }
+
+                    DscalVector(e, lp1, 1.0 / e[l]);
+                    e[lp1] = 1.0 + e[lp1];
+                }
+
+                e[l] = -e[l];
+                if (lp1 < matrixCopy.RowCount && e[l] != 0.0)
+                {
+                    // Apply the transformation.
+                    for (i = lp1; i < matrixCopy.RowCount; i++)
+                    {
+                        work[i] = 0.0;
+                    }
+
+                    for (j = lp1; j < matrixCopy.ColumnCount; j++)
+                    {
+                        for (var ii = lp1; ii < matrixCopy.RowCount; ii++)
+                        {
+                            work[ii] += e[j] * matrixCopy.At(ii, j);
+                        }
+                    }
+
+                    for (j = lp1; j < matrixCopy.ColumnCount; j++)
+                    {
+                        var ww = -e[j] / e[lp1];
+                        for (var ii = lp1; ii < matrixCopy.RowCount; ii++)
+                        {
+                            matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii]));
+                        }
+                    }
+                }
+
+                if (ComputeVectors)
+                {
+                    // Place the transformation in v for subsequent back multiplication.
+                    for (i = lp1; i < matrixCopy.ColumnCount; i++)
+                    {
+                        MatrixVT.At(i, l, e[i]);
+                    }
+                }
+            }
+
+            // Set up the final bidiagonal matrixCopy or order m.
+            var m = Math.Min(matrixCopy.ColumnCount, matrixCopy.RowCount + 1);
+            var nctp1 = nct + 1;
+            var nrtp1 = nrt + 1;
+            if (nct < matrixCopy.ColumnCount)
+            {
+                VectorS[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1));
+            }
+
+            if (matrixCopy.RowCount < m)
+            {
+                VectorS[m - 1] = 0.0;
+            }
+
+            if (nrtp1 < m)
+            {
+                e[nrtp1 - 1] = matrixCopy.At((nrtp1 - 1), (m - 1));
+            }
+
+            e[m - 1] = 0.0;
+
+            // If required, generate u.
+            if (ComputeVectors)
+            {
+                for (j = nctp1 - 1; j < ncu; j++)
+                {
+                    for (i = 0; i < matrixCopy.RowCount; i++)
+                    {
+                        MatrixU.At(i, j, 0.0);
+                    }
+
+                    MatrixU.At(j, j, 1.0);
+                }
+
+                for (l = nct - 1; l >= 0; l--)
+                {
+                    if (VectorS[l] != 0.0)
+                    {
+                        for (j = l + 1; j < ncu; j++)
+                        {
+                            t = -Ddot(MatrixU, matrixCopy.RowCount, l, j, l) / MatrixU.At(l, l);
+                            for (var ii = l; ii < matrixCopy.RowCount; ii++)
+                            {
+                                MatrixU.At(ii, j, MatrixU.At(ii, j) + (t * MatrixU.At(ii, l)));
+                            }
+                        }
+
+                        DscalColumn(MatrixU, matrixCopy.RowCount, l, l, -1.0);
+                        MatrixU.At(l, l, 1.0 + MatrixU.At(l, l));
+                        for (i = 0; i < l; i++)
+                        {
+                            MatrixU.At(i, l, 0.0);
+                        }
+                    }
+                    else
+                    {
+                        for (i = 0; i < matrixCopy.RowCount; i++)
+                        {
+                            MatrixU.At(i, l, 0.0);
+                        }
+
+                        MatrixU.At(l, l, 1.0);
+                    }
+                }
+            }
+
+            // If it is required, generate v.
+            if (ComputeVectors)
+            {
+                for (l = matrixCopy.ColumnCount - 1; l >= 0; l--)
+                {
+                    lp1 = l + 1;
+                    if (l < nrt)
+                    {
+                        if (e[l] != 0.0)
+                        {
+                            for (j = lp1; j < matrixCopy.ColumnCount; j++)
+                            {
+                                t = -Ddot(MatrixVT, matrixCopy.ColumnCount, l, j, lp1) / MatrixVT.At(lp1, l);
+                                for (var ii = l; ii < matrixCopy.ColumnCount; ii++)
+                                {
+                                    MatrixVT.At(ii, j, MatrixVT.At(ii, j) + (t * MatrixVT.At(ii, l)));
+                                }
+                            }
+                        }
+                    }
+
+                    for (i = 0; i < matrixCopy.ColumnCount; i++)
+                    {
+                        MatrixVT.At(i, l, 0.0);
+                    }
+
+                    MatrixVT.At(l, l, 1.0);
+                }
+            }
+
+            // Transform s and e so that they are  double .
+            for (i = 0; i < m; i++)
+            {
+                double r;
+                if (VectorS[i] != 0.0)
+                {
+                    t = VectorS[i];
+                    r = VectorS[i] / t;
+                    VectorS[i] = t;
+                    if (i < m - 1)
+                    {
+                        e[i] = e[i] / r;
+                    }
+
+                    if (ComputeVectors)
+                    {
+                        DscalColumn(MatrixU, matrixCopy.RowCount, i, 0, r);
+                    }
+                }
+
+                // Exit
+                if (i == m - 1)
+                {
+                    break;
+                }
+
+                if (e[i] != 0.0)
+                {
+                    t = e[i];
+                    r = t / e[i];
+                    e[i] = t;
+                    VectorS[i + 1] = VectorS[i + 1] * r;
+                    if (ComputeVectors)
+                    {
+                        DscalColumn(MatrixVT, matrixCopy.ColumnCount, i + 1, 0, r);
+                    }
+                }
+            }
+
+            // Main iteration loop for the singular values.
+            var mn = m;
+            var iter = 0;
+
+            while (m > 0)
+            {
+                // Quit if all the singular values have been found. If too many iterations have been performed, 
+                // throw exception that Convergence Failed
+                if (iter >= Maxiter)
+                {
+                     throw new ArgumentException(Resources.ConvergenceFailed);
+                }
+
+                // This section of the program inspects for negligible elements in the s and e arrays. On
+                // completion the variables kase and l are set as follows.
+                // Kase = 1     if VectorS[m] and e[l-1] are negligible and l < m
+                // Kase = 2     if VectorS[l] is negligible and l < m
+                // Kase = 3     if e[l-1] is negligible, l < m, and VectorS[l, ..., VectorS[m] are not negligible (qr step).
+                // Лase = 4     if e[m-1] is negligible (convergence).
+                double ztest;
+                double test;
+                for (l = m - 2; l >= 0; l--)
+                {
+                    test = Math.Abs(VectorS[l]) + Math.Abs(VectorS[l + 1]);
+                    ztest = test + Math.Abs(e[l]);
+                    if (ztest.AlmostEqualInDecimalPlaces(test, 15))
+                    {
+                        e[l] = 0.0;
+                        break;
+                    }
+                }
+
+                int kase;
+                if (l == m - 2)
+                {
+                    kase = 4;
+                }
+                else
+                {
+                    int ls;
+                    for (ls = m - 1; ls > l; ls--)
+                    {
+                        test = 0.0;
+                        if (ls != m - 1)
+                        {
+                            test = test + Math.Abs(e[ls]);
+                        }
+
+                        if (ls != l + 1)
+                        {
+                            test = test + Math.Abs(e[ls - 1]);
+                        }
+
+                        ztest = test + Math.Abs(VectorS[ls]);
+                        if (ztest.AlmostEqualInDecimalPlaces(test, 15))
+                        {
+                            VectorS[ls] = 0.0;
+                            break;
+                        }
+                    }
+
+                    if (ls == l)
+                    {
+                        kase = 3;
+                    }
+                    else if (ls == m - 1)
+                    {
+                        kase = 1;
+                    }
+                    else
+                    {
+                        kase = 2;
+                        l = ls;
+                    }
+                }
+
+                l = l + 1;
+
+                // Perform the task indicated by kase.
+                int k;
+                double f;
+                switch (kase)
+                {
+                    // Deflate negligible VectorS[m].
+                    case 1:
+                        f = e[m - 2];
+                        e[m - 2] = 0.0;
+                        double t1;
+                        for (var kk = l; kk < m - 1; kk++)
+                        {
+                            k = m - 2 - kk + l;
+                            t1 = VectorS[k];
+                            Drotg(ref t1, ref f, ref cs, ref sn);
+                            VectorS[k] = t1;
+                            if (k != l)
+                            {
+                                f = -sn * e[k - 1];
+                                e[k - 1] = cs * e[k - 1];
+                            }
+
+                            if (ComputeVectors)
+                            {
+                                Drot(MatrixVT, matrixCopy.ColumnCount, k, m - 1, cs, sn);
+                            }
+                        }
+
+                        break;
+
+                    // Split at negligible VectorS[l].
+                    case 2:
+                        f = e[l - 1];
+                        e[l - 1] = 0.0;
+                        for (k = l; k < m; k++)
+                        {
+                            t1 = VectorS[k];
+                            Drotg(ref t1, ref f, ref cs, ref sn);
+                            VectorS[k] = t1;
+                            f = -sn * e[k];
+                            e[k] = cs * e[k];
+                            if (ComputeVectors)
+                            {
+                                Drot(MatrixU, matrixCopy.RowCount, k, l - 1, cs, sn);
+                            }
+                        }
+
+                        break;
+
+                    // Perform one qr step.
+                    case 3:
+                        // Calculate the shift.
+                        var scale = 0.0;
+                        scale = Math.Max(scale, Math.Abs(VectorS[m - 1]));
+                        scale = Math.Max(scale, Math.Abs(VectorS[m - 2]));
+                        scale = Math.Max(scale, Math.Abs(e[m - 2]));
+                        scale = Math.Max(scale, Math.Abs(VectorS[l]));
+                        scale = Math.Max(scale, Math.Abs(e[l]));
+                        var sm = VectorS[m - 1] / scale;
+                        var smm1 = VectorS[m - 2] / scale;
+                        var emm1 = e[m - 2] / scale;
+                        var sl = VectorS[l] / scale;
+                        var el = e[l] / scale;
+                        var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0;
+                        var c = (sm * emm1) * (sm * emm1);
+                        var shift = 0.0;
+                        if (b != 0.0 || c != 0.0)
+                        {
+                            shift = Math.Sqrt((b * b) + c);
+                            if (b < 0.0)
+                            {
+                                shift = -shift;
+                            }
+
+                            shift = c / (b + shift);
+                        }
+
+                        f = ((sl + sm) * (sl - sm)) + shift;
+                        var g = sl * el;
+
+                        // Chase zeros.
+                        for (k = l; k < m - 1; k++)
+                        {
+                            Drotg(ref f, ref g, ref cs, ref sn);
+                            if (k != l)
+                            {
+                                e[k - 1] = f;
+                            }
+
+                            f = (cs * VectorS[k]) + (sn * e[k]);
+                            e[k] = (cs * e[k]) - (sn * VectorS[k]);
+                            g = sn * VectorS[k + 1];
+                            VectorS[k + 1] = cs * VectorS[k + 1];
+                            if (ComputeVectors)
+                            {
+                                Drot(MatrixVT, matrixCopy.ColumnCount, k, k + 1, cs, sn);
+                            }
+
+                            Drotg(ref f, ref g, ref cs, ref sn);
+                            VectorS[k] = f;
+                            f = (cs * e[k]) + (sn * VectorS[k + 1]);
+                            VectorS[k + 1] = (-sn * e[k]) + (cs * VectorS[k + 1]);
+                            g = sn * e[k + 1];
+                            e[k + 1] = cs * e[k + 1];
+                            if (ComputeVectors && k < matrixCopy.RowCount)
+                            {
+                                Drot(MatrixU, matrixCopy.RowCount, k, k + 1, cs, sn);
+                            }
+                        }
+
+                        e[m - 2] = f;
+                        iter = iter + 1;
+                        break;
+
+                    // Convergence.
+                    case 4:
+                        // Make the singular value  positive
+                        if (VectorS[l] < 0.0)
+                        {
+                            VectorS[l] = -VectorS[l];
+                            if (ComputeVectors)
+                            {
+                                DscalColumn(MatrixVT, matrixCopy.ColumnCount, l, 0, -1.0);
+                            }
+                        }
+
+                        // Order the singular value.
+                        while (l != mn - 1)
+                        {
+                            if (VectorS[l] >= VectorS[l + 1])
+                            {
+                                break;
+                            }
+
+                            t = VectorS[l];
+                            VectorS[l] = VectorS[l + 1];
+                            VectorS[l + 1] = t;
+                            if (ComputeVectors && l < matrixCopy.ColumnCount)
+                            {
+                                Dswap(MatrixVT, matrixCopy.ColumnCount, l, l + 1);
+                            }
+
+                            if (ComputeVectors && l < matrixCopy.RowCount)
+                            {
+                                Dswap(MatrixU, matrixCopy.RowCount, l, l + 1);
+                            }
+
+                            l = l + 1;
+                        }
+
+                        iter = 0;
+                        m = m - 1;
+                        break;
+                }
+            }
+
+            if (ComputeVectors)
+            {
+                MatrixVT = MatrixVT.Transpose();
+            }
+
+            // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses
+            // a singular vector of length mRows+1 when mRows < mColumns. The last element is not used and needs to be removed.
+            // we should port lapack's svd routine to remove this problem.
+            if (matrixCopy.RowCount < matrixCopy.ColumnCount)
+            {
+                nm--;
+                var tmp = matrixCopy.CreateVector(nm);
+                for (i = 0; i < nm; i++)
+                {
+                    tmp[i] = VectorS[i];
+                }
+
+                VectorS = tmp;
+            }
+        }
+
+        /// <summary>
+        /// Calculates absolute value of <paramref name="z1"/> multiplied on signum function of <paramref name="z2"/>
+        /// </summary>
+        /// <param name="z1">Double value z1</param>
+        /// <param name="z2">Double value z2</param>
+        /// <returns>Result multiplication of signum function and absolute value</returns>
+        private static double Dsign(double z1, double z2)
+        {
+            return Math.Abs(z1) * (z2 / Math.Abs(z2));
+        }
+
+        /// <summary>
+        /// Swap column  <paramref name="columnA"/>  and  <paramref name="columnB"/>
+        /// </summary>
+        /// <param name="a">Source matrix</param>
+        /// <param name="rowCount">The number of rows in <paramref name="a"/></param>
+        /// <param name="columnA">Column A index to swap</param>
+        /// <param name="columnB">Column B index to swap</param>
+        private static void Dswap(Matrix a, int rowCount, int columnA, int columnB)
+        {
+            for (var i = 0; i < rowCount; i++)
+            {
+                var z = a.At(i, columnA);
+                a.At(i, columnA, a.At(i, columnB));
+                a.At(i, columnB, z);
+            }
+        }
+
+        /// <summary>
+        /// Scale column <paramref name="column"/> by <paramref name="z"/> starting from row <paramref name="rowStart"/>
+        /// </summary>
+        /// <param name="a">Source matrix</param>
+        /// <param name="rowCount">The number of rows in <paramref name="a"/> </param>
+        /// <param name="column">Column to scale</param>
+        /// <param name="rowStart">Row to scale from</param>
+        /// <param name="z">Scale value</param>
+        private static void DscalColumn(Matrix a, int rowCount, int column, int rowStart, double z)
+        {
+            for (var i = rowStart; i < rowCount; i++)
+            {
+                a.At(i, column, a.At(i, column) * z);
+            }
+        }
+
+        /// <summary>
+        /// Scale vector <paramref name="a"/> by <paramref name="z"/> starting from index <paramref name="start"/>
+        /// </summary>
+        /// <param name="a">Source vector</param>
+        /// <param name="start">Row to scale from</param>
+        /// <param name="z">Scale value</param>
+        private static void DscalVector(double[] a, int start, double z)
+        {
+            for (var i = start; i < a.Length; i++)
+            {
+                a[i] = a[i] * z;
+            }
+        }
+
+        /// <summary>
+        /// Given the Cartesian coordinates (da, db) of a point p, these fucntion return the parameters da, db, c, and s 
+        /// associated with the Givens rotation that zeros the y-coordinate of the point.
+        /// </summary>
+        /// <param name="da">Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation</param>
+        /// <param name="db">Provides the y-coordinate of the point p. On exit contains the parameter z associated with the Givens rotation</param>
+        /// <param name="c">Contains the parameter c associated with the Givens rotation</param>
+        /// <param name="s">Contains the parameter s associated with the Givens rotation</param>
+        /// <remarks>This is equivalent to the DROTG LAPACK routine.</remarks>
+        private static void Drotg(ref double da, ref double db, ref double c, ref double s)
+        {
+            double r, z;
+
+            var roe = db;
+            var absda = Math.Abs(da);
+            var absdb = Math.Abs(db);
+            if (absda > absdb)
+            {
+                roe = da;
+            }
+
+            var scale = absda + absdb;
+            if (scale == 0.0)
+            {
+                c = 1.0;
+                s = 0.0;
+                r = 0.0;
+                z = 0.0;
+            }
+            else
+            {
+                var sda = da / scale;
+                var sdb = db / scale;
+                r = scale * Math.Sqrt((sda * sda) + (sdb * sdb));
+                if (roe < 0.0)
+                {
+                    r = -r;
+                }
+
+                c = da / r;
+                s = db / r;
+                z = 1.0;
+                if (absda > absdb)
+                {
+                    z = s;
+                }
+
+                if (absdb >= absda && c != 0.0)
+                {
+                    z = 1.0 / c;
+                }
+            }
+
+            da = r;
+            db = z;
+        }
+
+        /// <summary>dded
+        /// Calculate Norm 2 of the column <paramref name="column"/> in matrix <paramref name="a"/> starting from row <paramref name="rowStart"/>
+        /// </summary>
+        /// <param name="a">Source matrix</param>
+        /// <param name="rowCount">The number of rows in <paramref name="a"/></param>
+        /// <param name="column">Column index</param>
+        /// <param name="rowStart">Start row index</param>
+        /// <returns>Norm2 (Euclidean norm) of trhe column</returns>
+        private static double Dnrm2Column(Matrix a, int rowCount, int column, int rowStart)
+        {
+            double s = 0;
+            for (var i = rowStart; i < rowCount; i++)
+            {
+                s += a.At(i, column) * a.At(i, column);
+            }
+
+            return Math.Sqrt(s);
+        }
+
+        /// <summary>
+        /// Calculate Norm 2 of the vector <paramref name="a"/> starting from index <paramref name="rowStart"/>
+        /// </summary>
+        /// <param name="a">Source vector</param>
+        /// <param name="rowStart">Start index</param>
+        /// <returns>Norm2 (Euclidean norm) of the vector</returns>
+        private static double Dnrm2Vector(double[] a, int rowStart)
+        {
+            double s = 0;
+            for (var i = rowStart; i < a.Length; i++)
+            {
+                s += a[i] * a[i];
+            }
+
+            return Math.Sqrt(s);
+        }
+
+        /// <summary>
+        /// Calculate dot product of <paramref name="columnA"/> and <paramref name="columnB"/>
+        /// </summary>
+        /// <param name="a">Source matrix</param>
+        /// <param name="rowCount">The number of rows in <paramref name="a"/></param>
+        /// <param name="columnA">Index of column A</param>
+        /// <param name="columnB">Index of column B</param>
+        /// <param name="rowStart">Starting row index</param>
+        /// <returns>Dot product value</returns>
+        private static double Ddot(Matrix a, int rowCount, int columnA, int columnB, int rowStart)
+        {
+            var z = 0.0;
+            for (var i = rowStart; i < rowCount; i++)
+            {
+                z += a.At(i, columnB) * a.At(i, columnA);
+            }
+
+            return z;
+        }
+
+        /// <summary>
+        /// Performs rotation of points in the plane. Given two vectors x <paramref name="columnA"/> and y <paramref name="columnB"/>, 
+        /// each vector element of these vectors is replaced as follows: x(i) = c*x(i) + s*y(i); y(i) = c*y(i) - s*x(i)
+        /// </summary>
+        /// <param name="a">Source matrix</param>
+        /// <param name="rowCount">The number of rows in <paramref name="a"/></param>
+        /// <param name="columnA">Index of column A</param>
+        /// <param name="columnB">Index of column B</param>
+        /// <param name="c">Scalar "c" value</param>
+        /// <param name="s">Scalar "s" value</param>
+        private static void Drot(Matrix a, int rowCount, int columnA, int columnB, double c, double s)
+        {
+            for (var i = 0; i < rowCount; i++)
+            {
+                var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB));
+                var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA));
+                a.At(i, columnB, tmp);
+                a.At(i, columnA, z);
+            }
+        }
+ 
+        /// <summary>
+        /// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
+        /// </summary>
+        /// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
+        public override void Solve(Matrix input, Matrix result)
+        {
+            // Check for proper arguments.
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            if (!ComputeVectors)
+            {
+                throw new InvalidOperationException(Resources.SingularVectorsNotComputed);
+            }
+
+            // The solution X should have the same number of columns as B
+            if (input.ColumnCount != result.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+            }
+
+            // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
+            if (MatrixU.RowCount != input.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
+            }
+
+            // The solution X row dimension is equal to the column dimension of A
+            if (MatrixVT.ColumnCount != result.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
+            }
+
+            var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount);
+            var bn = input.ColumnCount;
+
+            var tmp = new double[MatrixVT.ColumnCount];
+
+            for (var k = 0; k < bn; k++)
+            {
+                for (var j = 0; j < MatrixVT.ColumnCount; j++)
+                {
+                    double value = 0;
+                    if (j < mn)
+                    {
+                        for (var i = 0; i < MatrixU.RowCount; i++)
+                        {
+                            value += MatrixU.At(i, j) * input.At(i, k);
+                        }
+
+                        value /= VectorS[j];
+                    }
+
+                    tmp[j] = value;
+                }
+
+                for (var j = 0; j < MatrixVT.ColumnCount; j++)
+                {
+                    double value = 0;
+                    for (var i = 0; i < MatrixVT.ColumnCount; i++)
+                    {
+                        value += MatrixVT.At(i, j) * tmp[i];
+                    }
+
+                    result[j, k] = value;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Solves a system of linear equations, <b>Ax = b</b>, with A SVD factorized.
+        /// </summary>
+        /// <param name="input">The right hand side vector, <b>b</b>.</param>
+        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>x</b>.</param>
+        public override void Solve(Vector input, Vector result)
+        {
+            if (input == null)
+            {
+                throw new ArgumentNullException("input");
+            }
+
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            if (!ComputeVectors)
+            {
+                throw new InvalidOperationException(Resources.SingularVectorsNotComputed);
+            }
+
+            // Ax=b where A is an m x n matrix
+            // Check that b is a column vector with m entries
+            if (MatrixU.RowCount != input.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            // Check that x is a column vector with n entries
+            if (MatrixVT.ColumnCount != result.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
+            }
+
+            var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount);
+            var tmp = new double[MatrixVT.ColumnCount];
+            double value;
+            for (var j = 0; j < MatrixVT.ColumnCount; j++)
+            {
+                value = 0;
+                if (j < mn)
+                {
+                    for (var i = 0; i < MatrixU.RowCount; i++)
+                    {
+                        value += MatrixU.At(i, j) * input[i];
+                    }
+
+                    value /= VectorS[j];
+                }
+
+                tmp[j] = value;
+            }
+
+            for (var j = 0; j < MatrixVT.ColumnCount; j++)
+            {
+                value = 0;
+                for (int i = 0; i < MatrixVT.ColumnCount; i++)
+                {
+                    value += MatrixVT.At(i, j) * tmp[i];
+                }
+
+                result[j] = value;
+            }
+        }
+    }
+}

src/Numerics/Numerics.csproj

@@ -103,6 +103,10 @@
     <Compile Include="Complex32.cs" />
     <Compile Include="LinearAlgebra\Double\Factorization\DenseQR.cs" />
     <Compile Include="LinearAlgebra\Double\Factorization\DenseSvd.cs" />
+    <Compile Include="LinearAlgebra\Double\Factorization\UserCholesky.cs" />
+    <Compile Include="LinearAlgebra\Double\Factorization\UserLU.cs" />
+    <Compile Include="LinearAlgebra\Double\Factorization\UserQR.cs" />
+    <Compile Include="LinearAlgebra\Double\Factorization\UserSvd.cs" />
     <Compile Include="LinearAlgebra\Double\ISolver.cs" />
     <Compile Include="LinearAlgebra\Double\Factorization\QR.cs" />
     <Compile Include="LinearAlgebra\Double\Factorization\Svd.cs" />

src/Silverlight/Silverlight.csproj

@@ -248,6 +248,18 @@
     <Compile Include="..\Numerics\LinearAlgebra\Double\Factorization\Svd.cs">
       <Link>LinearAlgebra\Double\Factorization\Svd.cs</Link>
     </Compile>
+    <Compile Include="..\Numerics\LinearAlgebra\Double\Factorization\UserCholesky.cs">
+      <Link>LinearAlgebra\Double\Factorization\UserCholesky.cs</Link>
+    </Compile>
+    <Compile Include="..\Numerics\LinearAlgebra\Double\Factorization\UserLU.cs">
+      <Link>LinearAlgebra\Double\Factorization\UserLU.cs</Link>
+    </Compile>
+    <Compile Include="..\Numerics\LinearAlgebra\Double\Factorization\UserQR.cs">
+      <Link>LinearAlgebra\Double\Factorization\UserQR.cs</Link>
+    </Compile>
+    <Compile Include="..\Numerics\LinearAlgebra\Double\Factorization\UserSvd.cs">
+      <Link>LinearAlgebra\Double\Factorization\UserSvd.cs</Link>
+    </Compile>
     <Compile Include="..\Numerics\LinearAlgebra\Double\ISolver.cs">
       <Link>LinearAlgebra\Double\ISolver.cs</Link>
     </Compile>

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

@@ -30,7 +30,6 @@
 
 namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
 {
-    using System.Collections.Generic;
     using MbUnit.Framework;
     using LinearAlgebra.Double;
     using LinearAlgebra.Double.Factorization;
@@ -44,23 +43,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         public void CanFactorizeIdentity(int order)
         {
             var I = DenseMatrix.Identity(order);
-            var C = I.Cholesky();
+            var factorC = I.Cholesky();
 
-            Assert.AreEqual(I.RowCount, C.Factor.RowCount);
-            Assert.AreEqual(I.ColumnCount, C.Factor.ColumnCount);
+            Assert.AreEqual(I.RowCount, factorC.Factor.RowCount);
+            Assert.AreEqual(I.ColumnCount, factorC.Factor.ColumnCount);
 
-            for (var i = 0; i < C.Factor.RowCount; i++)
+            for (var i = 0; i < factorC.Factor.RowCount; i++)
             {
-                for (var j = 0; j < C.Factor.ColumnCount; j++)
+                for (var j = 0; j < factorC.Factor.ColumnCount; j++)
                 {
-                    if (i == j)
-                    {
-                        Assert.AreEqual(1.0, C.Factor[i, j]);
-                    }
-                    else
-                    {
-                        Assert.AreEqual(0.0, C.Factor[i, j]);
-                    }
+                    Assert.AreEqual(i == j ? 1.0 : 0.0, factorC.Factor[i, j]);
                 }
             }
         }
@@ -71,7 +63,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         {
             var I = DenseMatrix.Identity(10);
             I[3, 3] = -4.0;
-            var C = I.Cholesky();
+            I.Cholesky();
         }
 
         [Test]
@@ -81,7 +73,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         public void CholeskyFailsWithNonSquareMatrix(int row, int col)
         {
             var I = new DenseMatrix(row, col);
-            var C = I.Cholesky();
+            I.Cholesky();
         }
 
         [Test]
@@ -91,9 +83,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         public void IdentityDeterminantIsOne(int order)
         {
             var I = DenseMatrix.Identity(order);
-            var C = I.Cholesky();
-            Assert.AreEqual(1.0, C.Determinant);
-            Assert.AreEqual(0.0, C.DeterminantLn);
+            var factorC = I.Cholesky();
+            Assert.AreEqual(1.0, factorC.Determinant);
+            Assert.AreEqual(0.0, factorC.DeterminantLn);
         }
 
         [Test]
@@ -106,30 +98,30 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanFactorizeRandomMatrix(int order)
         {
-            var X = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order);
-            var chol = X.Cholesky();
-            var C = chol.Factor;
+            var matrixX = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order);
+            var chol = matrixX.Cholesky();
+            var factorC = chol.Factor;
 
             // Make sure the Cholesky factor has the right dimensions.
-            Assert.AreEqual(order, C.RowCount);
-            Assert.AreEqual(order, C.ColumnCount);
+            Assert.AreEqual(order, factorC.RowCount);
+            Assert.AreEqual(order, factorC.ColumnCount);
 
             // Make sure the Cholesky factor is lower triangular.
-            for (int i = 0; i < C.RowCount; i++) 
+            for (var i = 0; i < factorC.RowCount; i++) 
             {
-                for (int j = i+1; j < C.ColumnCount; j++)
+                for (var j = i+1; j < factorC.ColumnCount; j++)
                 {
-                    Assert.AreEqual(0.0, C[i, j]);
+                    Assert.AreEqual(0.0, factorC[i, j]);
                 }
             }
 
             // Make sure the cholesky factor times it's transpose is the original matrix.
-            var XfromC = C * C.Transpose();
-            for (int i = 0; i < XfromC.RowCount; i++) 
+            var matrixXfromC = factorC * factorC.Transpose();
+            for (var i = 0; i < matrixXfromC.RowCount; i++) 
             {
-                for (int j = 0; j < XfromC.ColumnCount; j++)
+                for (var j = 0; j < matrixXfromC.ColumnCount; j++)
                 {
-                    Assert.AreApproximatelyEqual(X[i,j], XfromC[i, j], 1.0e-11);
+                    Assert.AreApproximatelyEqual(matrixX[i,j], matrixXfromC[i, j], 1.0e-11);
                 }
             }
         }
@@ -144,28 +136,28 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [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 matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var b = MatrixLoader.GenerateRandomDenseVector(order);
             var x = chol.Solve(b);
 
             Assert.AreEqual(b.Count, x.Count);
 
-            var bReconstruct = A * x;
+            var bReconstruct = matrixA * x;
 
             // Check the reconstruction.
-            for (int i = 0; i < order; i++)
+            for (var 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 (var i = 0; i < matrixA.RowCount; i++)
             {
-                for (int j = 0; j < A.ColumnCount; j++)
+                for (var j = 0; j < matrixA.ColumnCount; j++)
                 {
-                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                 }
             }
         }
@@ -180,32 +172,32 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [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);
+            var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(row);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, col);
+            var matrixX = chol.Solve(matrixB);
 
-            Assert.AreEqual(B.RowCount, X.RowCount);
-            Assert.AreEqual(B.ColumnCount, X.ColumnCount);
+            Assert.AreEqual(matrixB.RowCount, matrixX.RowCount);
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
 
-            var BReconstruct = A * X;
+            var matrixBReconstruct = matrixA * matrixX;
 
             // Check the reconstruction.
-            for (int i = 0; i < B.RowCount; i++)
+            for (var i = 0; i < matrixB.RowCount; i++)
             {
-                for (int j = 0; j < B.ColumnCount; j++)
+                for (var j = 0; j < matrixB.ColumnCount; j++)
                 {
-                    Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11);
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
                 }
             }
 
             // Make sure A didn't change.
-            for (int i = 0; i < A.RowCount; i++)
+            for (var i = 0; i < matrixA.RowCount; i++)
             {
-                for (int j = 0; j < A.ColumnCount; j++)
+                for (var j = 0; j < matrixA.ColumnCount; j++)
                 {
-                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                 }
             }
         }
@@ -220,35 +212,35 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [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 matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var b = MatrixLoader.GenerateRandomDenseVector(order);
             var bCopy = b.Clone();
             var x = new DenseVector(order);
             chol.Solve(b, x);
 
             Assert.AreEqual(b.Count, x.Count);
 
-            var bReconstruct = A * x;
+            var bReconstruct = matrixA * x;
 
             // Check the reconstruction.
-            for (int i = 0; i < order; i++)
+            for (var 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 (var i = 0; i < matrixA.RowCount; i++)
             {
-                for (int j = 0; j < A.ColumnCount; j++)
+                for (var j = 0; j < matrixA.ColumnCount; j++)
                 {
-                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                 }
             }
 
             // Make sure b didn't change.
-            for (int i = 0; i < order; i++)
+            for (var i = 0; i < order; i++)
             {
                 Assert.AreEqual(bCopy[i], b[i]);
             }
@@ -264,43 +256,43 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [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);
+            var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(row);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, col);
+            var matrixBCopy = matrixB.Clone();
+            var matrixX = new DenseMatrix(row, col);
+            chol.Solve(matrixB, matrixX);
 
-            Assert.AreEqual(B.RowCount, X.RowCount);
-            Assert.AreEqual(B.ColumnCount, X.ColumnCount);
+            Assert.AreEqual(matrixB.RowCount, matrixX.RowCount);
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
 
-            var BReconstruct = A * X;
+            var matrixBReconstruct = matrixA * matrixX;
 
             // Check the reconstruction.
-            for (int i = 0; i < B.RowCount; i++)
+            for (var i = 0; i < matrixB.RowCount; i++)
             {
-                for (int j = 0; j < B.ColumnCount; j++)
+                for (var j = 0; j < matrixB.ColumnCount; j++)
                 {
-                    Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11);
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
                 }
             }
 
             // Make sure A didn't change.
-            for (int i = 0; i < A.RowCount; i++)
+            for (var i = 0; i < matrixA.RowCount; i++)
             {
-                for (int j = 0; j < A.ColumnCount; j++)
+                for (var j = 0; j < matrixA.ColumnCount; j++)
                 {
-                    Assert.AreEqual(ACopy[i, j], A[i, j]);
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                 }
             }
 
             // Make sure B didn't change.
-            for (int i = 0; i < B.RowCount; i++)
+            for (var i = 0; i < matrixB.RowCount; i++)
             {
-                for (int j = 0; j < B.ColumnCount; j++)
+                for (var j = 0; j < matrixB.ColumnCount; j++)
                 {
-                    Assert.AreEqual(BCopy[i, j], B[i, j]);
+                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
                 }
             }
         }

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

@@ -101,7 +101,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanFactorizeRandomMatrix(int order)
         {
-            var matrixX = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixX = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var factorLU = matrixX.LU();
             var matrixL = factorLU.L;
             var matrixU = factorLU.U;
@@ -154,11 +154,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomVector(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorLU = matrixA.LU();
 
-            var vectorb = MatrixLoader.GenerateRandomVector(order);
+            var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
             var resultx = factorLU.Solve(vectorb);
 
             Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
@@ -191,11 +191,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomMatrix(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorLU = matrixA.LU();
 
-            var matrixB = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixX = factorLU.Solve(matrixB);
 
             // The solution X row dimension is equal to the column dimension of A
@@ -234,10 +234,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomVectorWhenResultVectorGiven(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorLU = matrixA.LU();
-            var vectorb = MatrixLoader.GenerateRandomVector(order);
+            var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
             var vectorbCopy = vectorb.Clone();
             var resultx = new DenseVector(order);
             factorLU.Solve(vectorb, resultx);
@@ -278,11 +278,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorLU = matrixA.LU();
 
-            var matrixB = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixBCopy = matrixB.Clone();
 
             var matrixX = new DenseMatrix(order, order);
@@ -333,7 +333,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanInverse(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorLU = matrixA.LU();
 

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

@@ -101,7 +101,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanFactorizeRandomMatrix(int row, int column)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(row, column);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column);
             var factorQR = matrixA.QR();
 
             // Make sure the R has the right dimensions.
@@ -145,11 +145,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomVector(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorQR = matrixA.QR();
 
-            var vectorb = MatrixLoader.GenerateRandomVector(order);
+            var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
             var resultx = factorQR.Solve(vectorb);
 
             Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
@@ -182,11 +182,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomMatrix(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorQR = matrixA.QR();
 
-            var matrixB = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixX = factorQR.Solve(matrixB);
 
             // The solution X row dimension is equal to the column dimension of A
@@ -225,10 +225,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomVectorWhenResultVectorGiven(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorQR = matrixA.QR();
-            var vectorb = MatrixLoader.GenerateRandomVector(order);
+            var vectorb = MatrixLoader.GenerateRandomDenseVector(order);
             var vectorbCopy = vectorb.Clone();
             var resultx = new DenseVector(order);
             factorQR.Solve(vectorb,resultx);
@@ -269,11 +269,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixACopy = matrixA.Clone();
             var factorQR = matrixA.QR();
 
-            var matrixB = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var matrixBCopy = matrixB.Clone();
 
             var matrixX = new DenseMatrix(order, order);

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

@@ -82,7 +82,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanFactorizeRandomMatrix(int row, int column)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(row, column);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column);
             var factorSvd = matrixA.Svd(true);
 
             // Make sure the U has the right dimensions.
@@ -115,7 +115,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CheckRankOfNonSquare(int row, int column)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(row, column);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column);
             var factorSvd = matrixA.Svd(true);
 
             var mn = Math.Min(row, column);
@@ -132,7 +132,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CheckRankSquare(int order)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(order, order);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order);
             var factorSvd = matrixA.Svd(true);
 
             if (factorSvd.Determinant != 0)
@@ -172,10 +172,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [ExpectedException(typeof(InvalidOperationException))]
         public void CannotSolveMatrixIfVectorsNotComputed()
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(10, 10);
             var factorSvd = matrixA.Svd(false);
 
-            var matrixB = MatrixLoader.GenerateRandomMatrix(10, 10);
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(10, 10);
             factorSvd.Solve(matrixB);
         }
 
@@ -183,10 +183,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [ExpectedException(typeof(InvalidOperationException))]
         public void CannotSolveVectorIfVectorsNotComputed()
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(10, 10);
             var factorSvd = matrixA.Svd(false);
 
-            var vectorb = MatrixLoader.GenerateRandomVector(10);
+            var vectorb = MatrixLoader.GenerateRandomDenseVector(10);
             factorSvd.Solve(vectorb);
         }
 
@@ -200,11 +200,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomVector(int row, int column)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(row, column);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column);
             var matrixACopy = matrixA.Clone();
             var factorSvd = matrixA.Svd(true);
 
-            var vectorb = MatrixLoader.GenerateRandomVector(row);
+            var vectorb = MatrixLoader.GenerateRandomDenseVector(row);
             var resultx = factorSvd.Solve(vectorb);
 
             Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
@@ -237,11 +237,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomMatrix(int row, int count)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(row, count);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, count);
             var matrixACopy = matrixA.Clone();
             var factorSvd = matrixA.Svd(true);
 
-            var matrixB = MatrixLoader.GenerateRandomMatrix(row, count);
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, count);
             var matrixX = factorSvd.Solve(matrixB);
 
             // The solution X row dimension is equal to the column dimension of A
@@ -280,10 +280,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomVectorWhenResultVectorGiven(int row, int column)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(row, column);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column);
             var matrixACopy = matrixA.Clone();
             var factorSvd = matrixA.Svd(true);
-            var vectorb = MatrixLoader.GenerateRandomVector(row);
+            var vectorb = MatrixLoader.GenerateRandomDenseVector(row);
             var vectorbCopy = vectorb.Clone();
             var resultx = new DenseVector(column);
             factorSvd.Solve(vectorb,resultx);
@@ -322,11 +322,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
         [MultipleAsserts]
         public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int column)
         {
-            var matrixA = MatrixLoader.GenerateRandomMatrix(row, column);
+            var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column);
             var matrixACopy = matrixA.Clone();
             var factorSvd = matrixA.Svd(true);
 
-            var matrixB = MatrixLoader.GenerateRandomMatrix(row, column);
+            var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, column);
             var matrixBCopy = matrixB.Clone();
 
             var matrixX = new DenseMatrix(column, column);

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

@@ -0,0 +1,299 @@
+// <copyright file="UserCholeskyTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebraTests.Double.Factorization
+{
+    using MbUnit.Framework;
+    using LinearAlgebra.Double.Factorization;
+
+    public class UserCholeskyTests
+    {
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        public void CanFactorizeIdentity(int order)
+        {
+            var I = UserDefinedMatrix.Identity(order);
+            var factorC = I.Cholesky();
+
+            Assert.AreEqual(I.RowCount, factorC.Factor.RowCount);
+            Assert.AreEqual(I.ColumnCount, factorC.Factor.ColumnCount);
+
+            for (var i = 0; i < factorC.Factor.RowCount; i++)
+            {
+                for (var j = 0; j < factorC.Factor.ColumnCount; j++)
+                {
+                    Assert.AreEqual(i == j ? 1.0 : 0.0, factorC.Factor[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [ExpectedArgumentException]
+        public void CholeskyFailsWithDiagonalNonPositiveDefiniteMatrix()
+        {
+            var I = UserDefinedMatrix.Identity(10);
+            I[3, 3] = -4.0;
+            I.Cholesky();
+        }
+
+        [Test]
+        [Row(3,5)]
+        [Row(5,3)]
+        [ExpectedArgumentException]
+        public void CholeskyFailsWithNonSquareMatrix(int row, int col)
+        {
+            var I = new UserDefinedMatrix(row, col);
+            I.Cholesky();
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        public void IdentityDeterminantIsOne(int order)
+        {
+            var I = UserDefinedMatrix.Identity(order);
+            var factorC = I.Cholesky();
+            Assert.AreEqual(1.0, factorC.Determinant);
+            Assert.AreEqual(0.0, factorC.DeterminantLn);
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanFactorizeRandomMatrix(int order)
+        {
+            var matrixX = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(order);
+            var chol = matrixX.Cholesky();
+            var factorC = chol.Factor;
+
+            // Make sure the Cholesky factor has the right dimensions.
+            Assert.AreEqual(order, factorC.RowCount);
+            Assert.AreEqual(order, factorC.ColumnCount);
+
+            // Make sure the Cholesky factor is lower triangular.
+            for (var i = 0; i < factorC.RowCount; i++) 
+            {
+                for (var j = i+1; j < factorC.ColumnCount; j++)
+                {
+                    Assert.AreEqual(0.0, factorC[i, j]);
+                }
+            }
+
+            // Make sure the cholesky factor times it's transpose is the original matrix.
+            var matrixXfromC = factorC * factorC.Transpose();
+            for (var i = 0; i < matrixXfromC.RowCount; i++) 
+            {
+                for (var j = 0; j < matrixXfromC.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixX[i,j], matrixXfromC[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 matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(order);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var b = MatrixLoader.GenerateRandomUserDefinedVector(order);
+            var x = chol.Solve(b);
+
+            Assert.AreEqual(b.Count, x.Count);
+
+            var bReconstruct = matrixA * x;
+
+            // Check the reconstruction.
+            for (var i = 0; i < order; i++)
+            {
+                Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[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 matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(row);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(row, col);
+            var matrixX = chol.Solve(matrixB);
+
+            Assert.AreEqual(matrixB.RowCount, matrixX.RowCount);
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomVectorWhenResultVectorGiven(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(order);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var b = MatrixLoader.GenerateRandomUserDefinedVector(order);
+            var bCopy = b.Clone();
+            var x = new UserDefinedVector(order);
+            chol.Solve(b, x);
+
+            Assert.AreEqual(b.Count, x.Count);
+
+            var bReconstruct = matrixA * x;
+
+            // Check the reconstruction.
+            for (var i = 0; i < order; i++)
+            {
+                Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure b didn't change.
+            for (var 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 matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(row);
+            var matrixACopy = matrixA.Clone();
+            var chol = matrixA.Cholesky();
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(row, col);
+            var matrixBCopy = matrixB.Clone();
+            var matrixX = new UserDefinedMatrix(row, col);
+            chol.Solve(matrixB, matrixX);
+
+            Assert.AreEqual(matrixB.RowCount, matrixX.RowCount);
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure B didn't change.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
+                }
+            }
+        }
+    }
+}

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

@@ -0,0 +1,363 @@
+// <copyright file="UserLUTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebraTests.Double.Factorization
+{
+    using MbUnit.Framework;
+    using LinearAlgebra.Double.Factorization;
+
+    public class UserLUTests
+    {
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        public void CanFactorizeIdentity(int order)
+        {
+            var matrixI = UserDefinedMatrix.Identity(order);
+            var factorLU = matrixI.LU();
+
+            // Check lower triangular part.
+            var matrixL = factorLU.L;
+            Assert.AreEqual(matrixI.RowCount, matrixL.RowCount);
+            Assert.AreEqual(matrixI.ColumnCount, matrixL.ColumnCount);
+            for (var i = 0; i < matrixL.RowCount; i++)
+            {
+                for (var j = 0; j < matrixL.ColumnCount; j++)
+                {
+                    Assert.AreEqual(i == j ? 1.0 : 0.0, matrixL[i, j]);
+                }
+            }
+
+            // Check upper triangular part.
+            var matrixU = factorLU.U;
+            Assert.AreEqual(matrixI.RowCount, matrixU.RowCount);
+            Assert.AreEqual(matrixI.ColumnCount, matrixU.ColumnCount);
+            for (var i = 0; i < matrixU.RowCount; i++)
+            {
+                for (var j = 0; j < matrixU.ColumnCount; j++)
+                {
+                    Assert.AreEqual(i == j ? 1.0 : 0.0, matrixU[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(3,5)]
+        [Row(5,3)]
+        [ExpectedArgumentException]
+        public void LUFailsWithNonSquareMatrix(int row, int col)
+        {
+            var I = new UserDefinedMatrix(row, col);
+            I.LU();
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        public void IdentityDeterminantIsOne(int order)
+        {
+            var I = UserDefinedMatrix.Identity(order);
+            var lu = I.LU();
+            Assert.AreEqual(1.0, lu.Determinant);
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanFactorizeRandomMatrix(int order)
+        {
+            var matrixX = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var factorLU = matrixX.LU();
+            var matrixL = factorLU.L;
+            var matrixU = factorLU.U;
+
+            // Make sure the factors have the right dimensions.
+            Assert.AreEqual(order, matrixL.RowCount);
+            Assert.AreEqual(order, matrixL.ColumnCount);
+            Assert.AreEqual(order, matrixU.RowCount);
+            Assert.AreEqual(order, matrixU.ColumnCount);
+
+            // Make sure the L factor is lower triangular.
+            for (var i = 0; i < matrixL.RowCount; i++) 
+            {
+                Assert.AreEqual(1.0, matrixL[i, i]);
+                for (var j = i+1; j < matrixL.ColumnCount; j++)
+                {
+                    Assert.AreEqual(0.0, matrixL[i, j]);
+                }
+            }
+
+            // Make sure the U factor is upper triangular.
+            for (var i = 0; i < matrixL.RowCount; i++)
+            {
+                for (var j = 0; j < i; j++)
+                {
+                    Assert.AreEqual(0.0, matrixU[i, j]);
+                }
+            }
+
+            // Make sure the LU factor times it's transpose is the original matrix.
+            var matrixXfromLU = matrixL * matrixU;
+            var permutationInverse = factorLU.P.Inverse();
+            matrixXfromLU.PermuteRows(permutationInverse);
+            for (var i = 0; i < matrixXfromLU.RowCount; i++) 
+            {
+                for (var j = 0; j < matrixXfromLU.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixX[i, j], matrixXfromLU[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 matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorLU = matrixA.LU();
+
+            var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order);
+            var resultx = factorLU.Solve(vectorb);
+
+            Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
+
+            var bReconstruct = matrixA * resultx;
+
+            // Check the reconstruction.
+            for (var i = 0; i < order; i++)
+            {
+                Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(4)]
+        [Row(8)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrix(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorLU = matrixA.LU();
+
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixX = factorLU.Solve(matrixB);
+
+            // The solution X row dimension is equal to the column dimension of A
+            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
+            // The solution X has the same number of columns as B
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomVectorWhenResultVectorGiven(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorLU = matrixA.LU();
+            var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order);
+            var vectorbCopy = vectorb.Clone();
+            var resultx = new UserDefinedVector(order);
+            factorLU.Solve(vectorb, resultx);
+
+            Assert.AreEqual(vectorb.Count, resultx.Count);
+
+            var bReconstruct = matrixA * resultx;
+
+            // Check the reconstruction.
+            for (var i = 0; i < vectorb.Count; i++)
+            {
+                Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure b didn't change.
+            for (var i = 0; i < vectorb.Count; i++)
+            {
+                Assert.AreEqual(vectorbCopy[i], vectorb[i]);
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(4)]
+        [Row(8)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorLU = matrixA.LU();
+
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixBCopy = matrixB.Clone();
+
+            var matrixX = new UserDefinedMatrix(order, order);
+            factorLU.Solve(matrixB, matrixX);
+
+            // The solution X row dimension is equal to the column dimension of A
+            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
+            // The solution X has the same number of columns as B
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure B didn't change.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(4)]
+        [Row(8)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanInverse(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorLU = matrixA.LU();
+
+            var matrixAInverse = factorLU.Inverse();
+
+            // The inverse dimension is equal A
+            Assert.AreEqual(matrixAInverse.RowCount, matrixAInverse.RowCount);
+            Assert.AreEqual(matrixAInverse.ColumnCount, matrixAInverse.ColumnCount);
+
+            var matrixIdentity = matrixA * matrixAInverse;
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Check if multiplication of A and AI produced identity matrix.
+            for (var i = 0; i < matrixIdentity.RowCount; i++)
+            {
+                Assert.AreApproximatelyEqual(matrixIdentity[i, i], 1.0, 1.0e-11);
+            }
+        }
+    }
+}

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

@@ -0,0 +1,316 @@
+// <copyright file="UserQRTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebraTests.Double.Factorization
+{
+    using MbUnit.Framework;
+    using LinearAlgebra.Double.Factorization;
+
+    public class UserQRTests
+    {
+
+        [Test]
+        [ExpectedArgumentNullException]
+        public void ConstructorNull()
+        {
+            new UserQR(null);
+        }
+
+        [Test]
+        [ExpectedArgumentException]
+        public void WideMatrixThrowsInvalidMatrixOperationException()
+        {
+            new UserQR(new UserDefinedMatrix(3, 4));
+        }
+        
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        public void CanFactorizeIdentity(int order)
+        {
+            var I = UserDefinedMatrix.Identity(order);
+            var factorQR = I.QR();
+
+            Assert.AreEqual(I.RowCount, factorQR.R.RowCount);
+            Assert.AreEqual(I.ColumnCount, factorQR.R.ColumnCount);
+
+            for (var i = 0; i < factorQR.R.RowCount; i++)
+            {
+                for (var j = 0; j < factorQR.R.ColumnCount; j++)
+                {
+                    if (i == j)
+                    {
+                        Assert.AreEqual(-1.0, factorQR.R[i, j]);
+                    }
+                    else
+                    {
+                        Assert.AreEqual(0.0, factorQR.R[i, j]);
+                    }
+                }
+            }
+        }
+
+      
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        public void IdentityDeterminantIsOne(int order)
+        {
+            var I = UserDefinedMatrix.Identity(order);
+            var factorQR = I.QR();
+            Assert.AreEqual(1.0, factorQR.Determinant);
+        }
+
+        [Test]
+        [Row(1,1)]
+        [Row(2,2)]
+        [Row(5,5)]
+        [Row(10,6)]
+        [Row(50,48)]
+        [Row(100,98)]
+        [MultipleAsserts]
+        public void CanFactorizeRandomMatrix(int row, int column)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column);
+            var factorQR = matrixA.QR();
+
+            // Make sure the R has the right dimensions.
+            Assert.AreEqual(row, factorQR.R.RowCount);
+            Assert.AreEqual(column, factorQR.R.ColumnCount);
+
+            // Make sure the Q has the right dimensions.
+            Assert.AreEqual(row, factorQR.Q.RowCount);
+            Assert.AreEqual(row, factorQR.Q.ColumnCount);
+
+            // Make sure the R factor is upper triangular.
+            for (var i = 0; i < factorQR.R.RowCount; i++) 
+            {
+                for (var j = 0; j < factorQR.R.ColumnCount; j++)
+                {
+                    if (i > j)
+                    {
+                        Assert.AreEqual(0.0, factorQR.R[i, j]);
+                    }
+                }
+            }
+
+            // Make sure the Q*R is the original matrix.
+            var matrixQfromR = factorQR.Q * factorQR.R;
+            for (int i = 0; i < matrixQfromR.RowCount; i++) 
+            {
+                for (int j = 0; j < matrixQfromR.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixA[i, j], matrixQfromR[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 matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorQR = matrixA.QR();
+
+            var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order);
+            var resultx = factorQR.Solve(vectorb);
+
+            Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
+
+            var bReconstruct = matrixA * resultx;
+
+            // Check the reconstruction.
+            for (var i = 0; i < order; i++)
+            {
+                Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(4)]
+        [Row(8)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrix(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorQR = matrixA.QR();
+
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixX = factorQR.Solve(matrixB);
+
+            // The solution X row dimension is equal to the column dimension of A
+            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
+            // The solution X has the same number of columns as B
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomVectorWhenResultVectorGiven(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorQR = matrixA.QR();
+            var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order);
+            var vectorbCopy = vectorb.Clone();
+            var resultx = new UserDefinedVector(order);
+            factorQR.Solve(vectorb,resultx);
+
+            Assert.AreEqual(vectorb.Count, resultx.Count);
+
+            var bReconstruct = matrixA * resultx;
+
+            // Check the reconstruction.
+            for (var i = 0; i < vectorb.Count; i++)
+            {
+                Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure b didn't change.
+            for (var i = 0; i < vectorb.Count; i++)
+            {
+                Assert.AreEqual(vectorbCopy[i], vectorb[i]);
+            }
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(4)]
+        [Row(8)]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixACopy = matrixA.Clone();
+            var factorQR = matrixA.QR();
+
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var matrixBCopy = matrixB.Clone();
+
+            var matrixX = new UserDefinedMatrix(order, order);
+            factorQR.Solve(matrixB,matrixX);
+
+            // The solution X row dimension is equal to the column dimension of A
+            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
+            // The solution X has the same number of columns as B
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure B didn't change.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
+                }
+            }
+        }
+    }
+}

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

@@ -0,0 +1,369 @@
+// <copyright file="UserSvdTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// 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
+// 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.LinearAlgebraTests.Double.Factorization
+{
+    using System;
+    using MbUnit.Framework;
+    using LinearAlgebra.Double.Factorization;
+
+    public class UserSvdTests
+    {
+
+        [Test]
+        [ExpectedArgumentNullException]
+        public void ConstructorNull()
+        {
+            new UserSvd(null, true);
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        public void CanFactorizeIdentity(int order)
+        {
+            var I = UserDefinedMatrix.Identity(order);
+            var factorSvd = I.Svd(true);
+
+            Assert.AreEqual(I.RowCount, factorSvd.U().RowCount);
+            Assert.AreEqual(I.RowCount, factorSvd.U().ColumnCount);
+
+            Assert.AreEqual(I.ColumnCount, factorSvd.VT().RowCount);
+            Assert.AreEqual(I.ColumnCount, factorSvd.VT().ColumnCount);
+
+            Assert.AreEqual(I.RowCount, factorSvd.W().RowCount);
+            Assert.AreEqual(I.ColumnCount, factorSvd.W().ColumnCount);
+
+            for (var i = 0; i < factorSvd.W().RowCount; i++)
+            {
+                for (var j = 0; j < factorSvd.W().ColumnCount; j++)
+                {
+                    Assert.AreEqual(i == j ? 1.0 : 0.0, factorSvd.W()[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1,1)]
+        [Row(2,2)]
+        [Row(5,5)]
+        [Row(10,6)]
+        [Row(48,52)]
+        [Row(100,93)]
+        [MultipleAsserts]
+        public void CanFactorizeRandomMatrix(int row, int column)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column);
+            var factorSvd = matrixA.Svd(true);
+
+            // Make sure the U has the right dimensions.
+            Assert.AreEqual(row, factorSvd.U().RowCount);
+            Assert.AreEqual(row, factorSvd.U().ColumnCount);
+
+            // Make sure the VT has the right dimensions.
+            Assert.AreEqual(column, factorSvd.VT().RowCount);
+            Assert.AreEqual(column, factorSvd.VT().ColumnCount);
+
+            // Make sure the W has the right dimensions.
+            Assert.AreEqual(row, factorSvd.W().RowCount);
+            Assert.AreEqual(column, factorSvd.W().ColumnCount);
+
+            // Make sure the U*W*VT is the original matrix.
+            var matrix = factorSvd.U() * factorSvd.W() * factorSvd.VT();
+            for (var i = 0; i < matrix.RowCount; i++) 
+            {
+                for (var j = 0; j < matrix.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixA[i, j], matrix[i, j], 1.0e-11);
+                }
+            }
+        }
+
+        [Test]
+        [Row(10, 8)]
+        [Row(48, 52)]
+        [Row(100, 93)]
+        [MultipleAsserts]
+        public void CheckRankOfNonSquare(int row, int column)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column);
+            var factorSvd = matrixA.Svd(true);
+
+            var mn = Math.Min(row, column);
+            Assert.AreEqual(factorSvd.Rank, mn);
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(5)]
+        [Row(9)]
+        [Row(50)]
+        [Row(90)]
+        [MultipleAsserts]
+        public void CheckRankSquare(int order)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order);
+            var factorSvd = matrixA.Svd(true);
+
+            if (factorSvd.Determinant != 0)
+            {
+                Assert.AreEqual(factorSvd.Rank, order);
+            }
+            else
+            {
+                Assert.AreEqual(factorSvd.Rank, order - 1);
+            }
+        }
+
+        [Test]
+        [Row(10)]
+        [Row(50)]
+        [Row(100)]
+        [MultipleAsserts]
+        public void CheckRankOfSquareSingular(int order)
+        {
+            var matrixA = new UserDefinedMatrix(order, order);
+            matrixA[0, 0] = 1;
+            matrixA[order - 1, order - 1] = 1;
+            for (var i = 1; i < order - 1; i++)
+            {
+                matrixA[i, i - 1] = 1;
+                matrixA[i, i + 1] = 1;
+                matrixA[i - 1, i] = 1;
+                matrixA[i + 1, i] = 1;
+            }
+            var factorSvd = matrixA.Svd(true);
+
+            Assert.AreEqual(factorSvd.Determinant, 0);
+            Assert.AreEqual(factorSvd.Rank, order - 1);
+        }
+
+        [Test]
+        [ExpectedException(typeof(InvalidOperationException))]
+        public void CannotSolveMatrixIfVectorsNotComputed()
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(10, 10);
+            var factorSvd = matrixA.Svd(false);
+
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(10, 10);
+            factorSvd.Solve(matrixB);
+        }
+
+        [Test]
+        [ExpectedException(typeof(InvalidOperationException))]
+        public void CannotSolveVectorIfVectorsNotComputed()
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(10, 10);
+            var factorSvd = matrixA.Svd(false);
+
+            var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(10);
+            factorSvd.Solve(vectorb);
+        }
+
+        [Test]
+        [Row(1, 1)]
+        [Row(2, 2)]
+        [Row(5, 5)]
+        [Row(9, 10)]
+        [Row(50, 50)]
+        [Row(90, 100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomVector(int row, int column)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column);
+            var matrixACopy = matrixA.Clone();
+            var factorSvd = matrixA.Svd(true);
+
+            var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(row);
+            var resultx = factorSvd.Solve(vectorb);
+
+            Assert.AreEqual(matrixA.ColumnCount, resultx.Count);
+
+            var bReconstruct = matrixA * resultx;
+
+            // Check the reconstruction.
+            for (var i = 0; i < vectorb.Count; i++)
+            {
+                Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1, 1)]
+        [Row(4, 4)]
+        [Row(7, 8)]
+        [Row(10, 10)]
+        [Row(45, 50)]
+        [Row(80, 100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrix(int row, int count)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, count);
+            var matrixACopy = matrixA.Clone();
+            var factorSvd = matrixA.Svd(true);
+
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(row, count);
+            var matrixX = factorSvd.Solve(matrixB);
+
+            // The solution X row dimension is equal to the column dimension of A
+            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
+            // The solution X has the same number of columns as B
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [Row(1, 1)]
+        [Row(2, 2)]
+        [Row(5, 5)]
+        [Row(9, 10)]
+        [Row(50, 50)]
+        [Row(90, 100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomVectorWhenResultVectorGiven(int row, int column)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column);
+            var matrixACopy = matrixA.Clone();
+            var factorSvd = matrixA.Svd(true);
+            var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(row);
+            var vectorbCopy = vectorb.Clone();
+            var resultx = new UserDefinedVector(column);
+            factorSvd.Solve(vectorb,resultx);
+
+            var bReconstruct = matrixA * resultx;
+
+            // Check the reconstruction.
+            for (var i = 0; i < vectorb.Count; i++)
+            {
+                Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11);
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure b didn't change.
+            for (var i = 0; i < vectorb.Count; i++)
+            {
+                Assert.AreEqual(vectorbCopy[i], vectorb[i]);
+            }
+        }
+
+        [Test]
+        [Row(1, 1)]
+        [Row(4, 4)]
+        [Row(7, 8)]
+        [Row(10, 10)]
+        [Row(45, 50)]
+        [Row(80, 100)]
+        [MultipleAsserts]
+        public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int column)
+        {
+            var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column);
+            var matrixACopy = matrixA.Clone();
+            var factorSvd = matrixA.Svd(true);
+
+            var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column);
+            var matrixBCopy = matrixB.Clone();
+
+            var matrixX = new UserDefinedMatrix(column, column);
+            factorSvd.Solve(matrixB,matrixX);
+
+            // The solution X row dimension is equal to the column dimension of A
+            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);
+            // The solution X has the same number of columns as B
+            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);
+
+            var matrixBReconstruct = matrixA * matrixX;
+
+            // Check the reconstruction.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11);
+                }
+            }
+
+            // Make sure A didn't change.
+            for (var i = 0; i < matrixA.RowCount; i++)
+            {
+                for (var j = 0; j < matrixA.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
+                }
+            }
+
+            // Make sure B didn't change.
+            for (var i = 0; i < matrixB.RowCount; i++)
+            {
+                for (var j = 0; j < matrixB.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]);
+                }
+            }
+        }
+    }
+}

src/UnitTests/LinearAlgebraTests/Double/MatrixLoader.cs

@@ -63,7 +63,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             }
         }
 
-        public static Matrix GenerateRandomMatrix(int row, int col)
+        public static Matrix GenerateRandomDenseMatrix(int row, int col)
         {
             // Fill a matrix with standard random numbers.
             var normal = new Distributions.Normal();
@@ -81,7 +81,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             return A;
         }
 
-        public static Matrix GenerateRandomPositiveDefiniteMatrix(int order)
+        public static Matrix GenerateRandomPositiveDefiniteDenseMatrix(int order)
         {
             // Fill a matrix with standard random numbers.
             var normal = new Distributions.Normal();
@@ -99,7 +99,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             return A.Transpose() * A;
         }
 
-        public static Vector GenerateRandomVector(int order)
+        public static Vector GenerateRandomDenseVector(int order)
         {
             // Fill a matrix with standard random numbers.
             var normal = new Distributions.Normal();
@@ -113,5 +113,56 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             // Generate a matrix which is positive definite.
             return v;
         }
+
+        public static Matrix GenerateRandomUserDefinedMatrix(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 UserDefinedMatrix(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 GenerateRandomPositiveDefiniteUserDefinedMatrix(int order)
+        {
+            // Fill a matrix with standard random numbers.
+            var normal = new Distributions.Normal();
+            normal.RandomSource = new Random.MersenneTwister(1);
+            var A = new UserDefinedMatrix(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 GenerateRandomUserDefinedVector(int order)
+        {
+            // Fill a matrix with standard random numbers.
+            var normal = new Distributions.Normal();
+            normal.RandomSource = new Random.MersenneTwister(1);
+            var v = new UserDefinedVector(order);
+            for (int i = 0; i < order; i++)
+            {
+                v[i] = normal.Sample();
+            }
+
+            // Generate a matrix which is positive definite.
+            return v;
+        }
     }
 }

src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrixTests.cs

@@ -36,6 +36,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
     {
         private readonly double[,] _data;
 
+        public UserDefinedMatrix(int order): base(order, order)
+        {
+            _data = new double[order, order];
+        }
+
         public UserDefinedMatrix(int rows, int columns) : base(rows, columns)
         {
             _data = new double[rows, columns];
@@ -65,6 +70,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         {
             return new UserDefinedVector(size);
         }
+
+        public static UserDefinedMatrix Identity(int order)
+        {
+            var m = new UserDefinedMatrix(order, order);
+            for (var i = 0; i < order; i++)
+            {
+                m[i, i] = 1.0;
+            }
+
+            return m;
+        }
     }
 
     public class UserDefinedMatrixTests : MatrixTests

src/UnitTests/UnitTests.csproj

@@ -90,6 +90,10 @@
     <Compile Include="ComplexTests\ComplexTest.cs" />
     <Compile Include="ComplexTests\Complex32Test.TextHandling.cs" />
     <Compile Include="ComplexTests\Complex32Test.cs" />
+    <Compile Include="LinearAlgebraTests\Double\Factorization\UserCholeskyTests.cs" />
+    <Compile Include="LinearAlgebraTests\Double\Factorization\UserLUTests.cs" />
+    <Compile Include="LinearAlgebraTests\Double\Factorization\UserQRTests.cs" />
+    <Compile Include="LinearAlgebraTests\Double\Factorization\UserSvdTests.cs" />
     <Compile Include="LinearAlgebraTests\Double\Factorization\SvdTests.cs" />
     <Compile Include="LinearAlgebraTests\Double\Factorization\QRTests.cs" />
     <Compile Include="LinearAlgebraTests\Double\SparseMatrixTests.cs" />