diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index 26314d10..4e3b52eb 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -1016,7 +1016,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
         public override LU<Complex> LU()
         {
-            return new DenseLU(this);
+            return DenseLU.Create(this);
         }
 
         public override QR<Complex> QR(QRMethod method = QRMethod.Thin)
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
index dc39c99e..b9a6e747 100644
--- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 {
@@ -38,6 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     using Numerics;
 #else
     using System.Numerics;
+
 #endif
 
     /// <summary>
@@ -48,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class DenseLU : LU
+    public sealed class DenseLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseLU"/> class. This object will compute the
@@ -57,7 +58,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// <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 DenseLU(DenseMatrix matrix)
+        public static DenseLU Create(DenseMatrix matrix)
         {
             if (matrix == null)
             {
@@ -70,12 +71,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             }
 
             // Create an array for the pivot indices.
-            Pivots = new int[matrix.RowCount];
+            var pivots = new int[matrix.RowCount];
 
             // Create a new matrix for the LU factors, then perform factorization (while overwriting).
-            var factors = (DenseMatrix)matrix.Clone();
-            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, Pivots);
-            Factors = factors;
+            var factors = (DenseMatrix) matrix.Clone();
+            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, pivots);
+
+            return new DenseLU(factors, pivots);
+        }
+
+        DenseLU(Matrix<Complex> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -128,7 +135,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
+            var dfactors = (DenseMatrix) Factors;
             Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
@@ -177,7 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
+            var dfactors = (DenseMatrix) Factors;
             Control.LinearAlgebraProvider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
@@ -187,7 +194,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// <returns>The inverse of this matrix.</returns>
         public override Matrix<Complex> Inverse()
         {
-            var result = (DenseMatrix)Factors.Clone();
+            var result = (DenseMatrix) Factors.Clone();
             Control.LinearAlgebraProvider.LUInverseFactored(result.Values, result.RowCount, Pivots);
             return result;
         }
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs
index 8d98abff..5e117c1b 100644
--- a/src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/LU.cs
@@ -3,7 +3,9 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+//
+// Copyright (c) 2009-2013 Math.NET
+//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,8 +14,10 @@
 // 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
@@ -46,6 +50,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     /// </remarks>
     public abstract class LU : LU<Complex>
     {
+        protected LU(Matrix<Complex> factors, int[] pivots)
+            : base(factors, pivots)
+        {
+        }
+
         /// <summary>
         /// Gets the determinant of the matrix for which the LU factorization was computed.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs
index 60ff3ca7..fd3f5763 100644
--- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 {
@@ -38,6 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     using Numerics;
 #else
     using System.Numerics;
+
 #endif
 
     /// <summary>
@@ -48,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class UserLU : LU
+    public sealed class UserLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="UserLU"/> class. This object will compute the
@@ -57,7 +58,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// <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<Complex> matrix)
+        public static UserLU Create(Matrix<Complex> matrix)
         {
             if (matrix == null)
             {
@@ -71,13 +72,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 
             // Create an array for the pivot indices.
             var order = matrix.RowCount;
-            Factors = matrix.Clone();
-            Pivots = new int[order];
-            
+            var factors = matrix.Clone();
+            var pivots = new int[order];
+
             // Initialize the pivot matrix to the identity permutation.
             for (var i = 0; i < order; i++)
             {
-                Pivots[i] = i;
+                pivots[i] = i;
             }
 
             var vectorLUcolj = new Complex[order];
@@ -86,7 +87,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 // Make a copy of the j-th column to localize references.
                 for (var i = 0; i < order; i++)
                 {
-                    vectorLUcolj[i] = Factors.At(i, j);
+                    vectorLUcolj[i] = factors.At(i, j);
                 }
 
                 // Apply previous transformations.
@@ -96,11 +97,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                     var s = Complex.Zero;
                     for (var k = 0; k < kmax; k++)
                     {
-                        s += Factors.At(i, k) * vectorLUcolj[k];
+                        s += factors.At(i, k)*vectorLUcolj[k];
                     }
 
                     vectorLUcolj[i] -= s;
-                    Factors.At(i, j, vectorLUcolj[i]);
+                    factors.At(i, j, vectorLUcolj[i]);
                 }
 
                 // Find pivot and exchange if necessary.
@@ -117,23 +118,30 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 {
                     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);
+                        var temp = factors.At(p, k);
+                        factors.At(p, k, factors.At(j, k));
+                        factors.At(j, k, temp);
                     }
 
-                    Pivots[j] = p;
+                    pivots[j] = p;
                 }
 
                 // Compute multipliers.
-                if (j < order & Factors.At(j, j) != 0.0)
+                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)));
+                        factors.At(i, j, (factors.At(i, j)/factors.At(j, j)));
                     }
                 }
             }
+
+            return new UserLU(factors, pivots);
+        }
+
+        UserLU(Matrix<Complex> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -189,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             }
 
             var order = Factors.RowCount;
-            
+
             // Solve L*Y = P*B
             for (var k = 0; k < order; k++)
             {
@@ -197,7 +205,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 {
                     for (var j = 0; j < result.ColumnCount; j++)
                     {
-                        var temp = result.At(k, j) * Factors.At(i, k);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -208,14 +216,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             {
                 for (var j = 0; j < result.ColumnCount; j++)
                 {
-                    result.At(k, j, (result.At(k, j) / Factors.At(k, k)));
+                    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);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -265,7 +273,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 result[p] = result[i];
                 result[i] = temp;
             }
-            
+
             var order = Factors.RowCount;
 
             // Solve L*Y = P*B
@@ -273,7 +281,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             {
                 for (var i = k + 1; i < order; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
 
@@ -283,7 +291,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 result[k] /= Factors.At(k, k);
                 for (var i = 0; i < k; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
index 5a6231b6..6690c7ed 100644
--- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
@@ -465,7 +465,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
         public override LU<Complex> LU()
         {
-            return new UserLU(this);
+            return UserLU.Create(this);
         }
 
         public override QR<Complex> QR(QRMethod method = QRMethod.Thin)
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index d2fd2239..79dcb352 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -1011,7 +1011,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
         public override LU<Complex32> LU()
         {
-            return new DenseLU(this);
+            return DenseLU.Create(this);
         }
 
         public override QR<Complex32> QR(QRMethod method = QRMethod.Thin)
diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs
index 71c8cca8..847bd3ea 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 {
@@ -43,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class DenseLU : LU
+    public sealed class DenseLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseLU"/> class. This object will compute the
@@ -52,7 +52,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// <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 DenseLU(DenseMatrix matrix)
+        public static DenseLU Create(DenseMatrix matrix)
         {
             if (matrix == null)
             {
@@ -65,12 +65,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             }
 
             // Create an array for the pivot indices.
-            Pivots = new int[matrix.RowCount];
+            var pivots = new int[matrix.RowCount];
 
             // Create a new matrix for the LU factors, then perform factorization (while overwriting).
-            var factors = (DenseMatrix)matrix.Clone();
-            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, Pivots);
-            Factors = factors;
+            var factors = (DenseMatrix) matrix.Clone();
+            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, pivots);
+
+            return new DenseLU(factors, pivots);
+        }
+
+        DenseLU(Matrix<Complex32> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -123,7 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
+            var dfactors = (DenseMatrix) Factors;
             Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
@@ -172,7 +178,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
+            var dfactors = (DenseMatrix) Factors;
             Control.LinearAlgebraProvider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
@@ -182,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// <returns>The inverse of this matrix.</returns>
         public override Matrix<Complex32> Inverse()
         {
-            var result = (DenseMatrix)Factors.Clone();
+            var result = (DenseMatrix) Factors.Clone();
             Control.LinearAlgebraProvider.LUInverseFactored(result.Values, result.RowCount, Pivots);
             return result;
         }
diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs
index 573392e3..349d1216 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/LU.cs
@@ -3,7 +3,9 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+//
+// Copyright (c) 2009-2013 Math.NET
+//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,8 +14,10 @@
 // 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
@@ -42,6 +46,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
     /// </remarks>
     public abstract class LU : LU<Complex32>
     {
+        protected LU(Matrix<Complex32> factors, int[] pivots)
+            : base(factors, pivots)
+        {
+        }
+
         /// <summary>
         /// Gets the determinant of the matrix for which the LU factorization was computed.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs
index 2e09244c..d87cdf5c 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 {
@@ -43,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class UserLU : LU
+    public sealed class UserLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="UserLU"/> class. This object will compute the
@@ -52,7 +52,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// <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<Complex32> matrix)
+        public static UserLU Create(Matrix<Complex32> matrix)
         {
             if (matrix == null)
             {
@@ -66,13 +66,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 
             // Create an array for the pivot indices.
             var order = matrix.RowCount;
-            Factors = matrix.Clone();
-            Pivots = new int[order];
-            
+            var factors = matrix.Clone();
+            var pivots = new int[order];
+
             // Initialize the pivot matrix to the identity permutation.
             for (var i = 0; i < order; i++)
             {
-                Pivots[i] = i;
+                pivots[i] = i;
             }
 
             var vectorLUcolj = new Complex32[order];
@@ -81,7 +81,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 // Make a copy of the j-th column to localize references.
                 for (var i = 0; i < order; i++)
                 {
-                    vectorLUcolj[i] = Factors.At(i, j);
+                    vectorLUcolj[i] = factors.At(i, j);
                 }
 
                 // Apply previous transformations.
@@ -91,11 +91,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                     var s = Complex32.Zero;
                     for (var k = 0; k < kmax; k++)
                     {
-                        s += Factors.At(i, k) * vectorLUcolj[k];
+                        s += factors.At(i, k)*vectorLUcolj[k];
                     }
 
                     vectorLUcolj[i] -= s;
-                    Factors.At(i, j, vectorLUcolj[i]);
+                    factors.At(i, j, vectorLUcolj[i]);
                 }
 
                 // Find pivot and exchange if necessary.
@@ -112,23 +112,30 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 {
                     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);
+                        var temp = factors.At(p, k);
+                        factors.At(p, k, factors.At(j, k));
+                        factors.At(j, k, temp);
                     }
 
-                    Pivots[j] = p;
+                    pivots[j] = p;
                 }
 
                 // Compute multipliers.
-                if (j < order & Factors.At(j, j) != 0.0f)
+                if (j < order & factors.At(j, j) != 0.0f)
                 {
                     for (var i = j + 1; i < order; i++)
                     {
-                        Factors.At(i, j, (Factors.At(i, j) / Factors.At(j, j)));
+                        factors.At(i, j, (factors.At(i, j)/factors.At(j, j)));
                     }
                 }
             }
+
+            return new UserLU(factors, pivots);
+        }
+
+        UserLU(Matrix<Complex32> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -184,7 +191,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             }
 
             var order = Factors.RowCount;
-            
+
             // Solve L*Y = P*B
             for (var k = 0; k < order; k++)
             {
@@ -192,7 +199,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 {
                     for (var j = 0; j < result.ColumnCount; j++)
                     {
-                        var temp = result.At(k, j) * Factors.At(i, k);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -203,14 +210,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             {
                 for (var j = 0; j < result.ColumnCount; j++)
                 {
-                    result.At(k, j, (result.At(k, j) / Factors.At(k, k)));
+                    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);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -260,7 +267,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 result[p] = result[i];
                 result[i] = temp;
             }
-            
+
             var order = Factors.RowCount;
 
             // Solve L*Y = P*B
@@ -268,7 +275,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             {
                 for (var i = k + 1; i < order; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
 
@@ -278,7 +285,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 result[k] /= Factors.At(k, k);
                 for (var i = 0; i < k; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
index 9bf95c5c..c4123ee1 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
@@ -460,7 +460,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
         public override LU<Complex32> LU()
         {
-            return new UserLU(this);
+            return UserLU.Create(this);
         }
 
         public override QR<Complex32> QR(QRMethod method = QRMethod.Thin)
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index a874d9bc..635226d8 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -1042,7 +1042,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
         public override LU<double> LU()
         {
-            return new DenseLU(this);
+            return DenseLU.Create(this);
         }
 
         public override QR<double> QR(QRMethod method = QRMethod.Thin)
diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs
index 55b788e2..64f4c5e0 100644
--- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs
+++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 {
@@ -41,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class DenseLU : LU
+    public sealed class DenseLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseLU"/> class. This object will compute the
@@ -50,7 +50,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// <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 DenseLU(DenseMatrix matrix)
+        public static DenseLU Create(DenseMatrix matrix)
         {
             if (matrix == null)
             {
@@ -63,12 +63,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             }
 
             // Create an array for the pivot indices.
-            Pivots = new int[matrix.RowCount];
+            var pivots = new int[matrix.RowCount];
 
             // Create a new matrix for the LU factors, then perform factorization (while overwriting).
-            var factors = (DenseMatrix)matrix.Clone();
-            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, Pivots);
-            Factors = factors;
+            var factors = (DenseMatrix) matrix.Clone();
+            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, pivots);
+
+            return new DenseLU(factors, pivots);
+        }
+
+        DenseLU(Matrix<double> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -118,11 +124,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             }
 
             // Copy the contents of input to result.
-            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble);
+            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
-            Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots,  dresult.Values);
+            var dfactors = (DenseMatrix) Factors;
+            Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
         /// <summary>
@@ -167,10 +173,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             }
 
             // Copy the contents of input to result.
-            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble);
+            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
+            var dfactors = (DenseMatrix) Factors;
             Control.LinearAlgebraProvider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
@@ -180,7 +186,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// <returns>The inverse of this matrix.</returns>
         public override Matrix<double> Inverse()
         {
-            var result = (DenseMatrix)Factors.Clone();
+            var result = (DenseMatrix) Factors.Clone();
             Control.LinearAlgebraProvider.LUInverseFactored(result.Values, result.RowCount, Pivots);
             return result;
         }
diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs
index 4f9c09f2..0e5f6f5d 100644
--- a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs
+++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs
@@ -3,7 +3,9 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+//
+// Copyright (c) 2009-2013 Math.NET
+//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,8 +14,10 @@
 // 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
@@ -40,6 +44,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
     /// </remarks>
     public abstract class LU : LU<double>
     {
+        protected LU(Matrix<double> factors, int[] pivots)
+            : base(factors, pivots)
+        {
+        }
+
         /// <summary>
         /// Gets the determinant of the matrix for which the LU factorization was computed.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs
index 3eb39276..9d3ffa6f 100644
--- a/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs
+++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 {
@@ -41,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class UserLU : LU
+    public sealed class UserLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="UserLU"/> class. This object will compute the
@@ -50,7 +50,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// <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<double> matrix)
+        public static UserLU Create(Matrix<double> matrix)
         {
             if (matrix == null)
             {
@@ -64,13 +64,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 
             // Create an array for the pivot indices.
             var order = matrix.RowCount;
-            Factors = matrix.Clone();
-            Pivots = new int[order];
-            
+            var factors = matrix.Clone();
+            var pivots = new int[order];
+
             // Initialize the pivot matrix to the identity permutation.
             for (var i = 0; i < order; i++)
             {
-                Pivots[i] = i;
+                pivots[i] = i;
             }
 
             var vectorLUcolj = new double[order];
@@ -79,7 +79,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 // Make a copy of the j-th column to localize references.
                 for (var i = 0; i < order; i++)
                 {
-                    vectorLUcolj[i] = Factors.At(i, j);
+                    vectorLUcolj[i] = factors.At(i, j);
                 }
 
                 // Apply previous transformations.
@@ -89,11 +89,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                     var s = 0.0;
                     for (var k = 0; k < kmax; k++)
                     {
-                        s += Factors.At(i, k) * vectorLUcolj[k];
+                        s += factors.At(i, k)*vectorLUcolj[k];
                     }
 
                     vectorLUcolj[i] -= s;
-                    Factors.At(i, j, vectorLUcolj[i]);
+                    factors.At(i, j, vectorLUcolj[i]);
                 }
 
                 // Find pivot and exchange if necessary.
@@ -110,23 +110,30 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 {
                     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);
+                        var temp = factors.At(p, k);
+                        factors.At(p, k, factors.At(j, k));
+                        factors.At(j, k, temp);
                     }
 
-                    Pivots[j] = p;
+                    pivots[j] = p;
                 }
 
                 // Compute multipliers.
-                if (j < order & Factors.At(j, j) != 0.0)
+                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)));
+                        factors.At(i, j, (factors.At(i, j)/factors.At(j, j)));
                     }
                 }
             }
+
+            return new UserLU(factors, pivots);
+        }
+
+        UserLU(Matrix<double> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -182,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             }
 
             var order = Factors.RowCount;
-            
+
             // Solve L*Y = P*B
             for (var k = 0; k < order; k++)
             {
@@ -190,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 {
                     for (var j = 0; j < result.ColumnCount; j++)
                     {
-                        var temp = result.At(k, j) * Factors.At(i, k);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -201,14 +208,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             {
                 for (var j = 0; j < result.ColumnCount; j++)
                 {
-                    result.At(k, j, (result.At(k, j) / Factors.At(k, k)));
+                    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);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -258,7 +265,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 result[p] = result[i];
                 result[i] = temp;
             }
-            
+
             var order = Factors.RowCount;
 
             // Solve L*Y = P*B
@@ -266,7 +273,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             {
                 for (var i = k + 1; i < order; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
 
@@ -276,7 +283,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 result[k] /= Factors.At(k, k);
                 for (var i = 0; i < k; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs
index 51d0e295..7440e30c 100644
--- a/src/Numerics/LinearAlgebra/Double/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs
@@ -466,7 +466,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
         public override LU<double> LU()
         {
-            return new UserLU(this);
+            return UserLU.Create(this);
         }
 
         public override QR<double> QR(QRMethod method = QRMethod.Thin)
diff --git a/src/Numerics/LinearAlgebra/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Factorization/LU.cs
index db569eea..85e581b4 100644
--- a/src/Numerics/LinearAlgebra/Factorization/LU.cs
+++ b/src/Numerics/LinearAlgebra/Factorization/LU.cs
@@ -3,7 +3,9 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+//
+// Copyright (c) 2009-2013 Math.NET
+//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,8 +14,10 @@
 // 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
@@ -40,60 +44,59 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization
     /// </remarks>
     /// <typeparam name="T">Supported data types are double, single, <see cref="Complex"/>, and <see cref="Complex32"/>.</typeparam>
     public abstract class LU<T> : ISolver<T>
-    where T : struct, IEquatable<T>, IFormattable
+        where T : struct, IEquatable<T>, IFormattable
     {
-        /// <summary>
-        /// Value of one for T.
-        /// </summary>
-        private static readonly T One = Builder<T>.Instance.One;
+        static readonly T One = Builder<T>.Instance.One;
 
-        /// <summary>
-        /// Gets or sets both the L and U factors in the same matrix.
-        /// </summary>
-        protected Matrix<T> Factors { get; set; }
+        readonly Lazy<Matrix<T>> _lazyL;
+        readonly Lazy<Matrix<T>> _lazyU;
+        readonly Lazy<Permutation> _lazyP;
 
-        /// <summary>
-        /// Gets or sets the pivot indices of the LU factorization.
-        /// </summary>
-        protected int[] Pivots { get; set; }
+        protected readonly Matrix<T> Factors;
+        protected readonly int[] Pivots;
+
+        protected LU(Matrix<T> factors, int[] pivots)
+        {
+            Factors = factors;
+            Pivots = pivots;
+
+            _lazyL = new Lazy<Matrix<T>>(ComputeL);
+            _lazyU = new Lazy<Matrix<T>>(Factors.UpperTriangle);
+            _lazyP = new Lazy<Permutation>(() => Permutation.FromInversions(Pivots));
+        }
+
+        Matrix<T> ComputeL()
+        {
+            var result = Factors.LowerTriangle();
+            for (var i = 0; i < result.RowCount; i++)
+            {
+                result.At(i, i, One);
+            }
+            return result;
+        }
 
         /// <summary>
         /// Gets the lower triangular factor.
         /// </summary>
-        public virtual Matrix<T> L
+        public Matrix<T> L
         {
-            get
-            {
-                var result = Factors.LowerTriangle();
-                for (var i = 0; i < result.RowCount; i++)
-                {
-                    result.At(i, i, One);
-                }
-
-                return result;
-            }
+            get { return _lazyL.Value; }
         }
 
         /// <summary>
         /// Gets the upper triangular factor.
         /// </summary>
-        public virtual Matrix<T> U
+        public Matrix<T> U
         {
-            get
-            {
-                return Factors.UpperTriangle();
-            }
+            get { return _lazyU.Value; }
         }
 
         /// <summary>
         /// Gets the permutation applied to LU factorization.
         /// </summary>
-        public virtual Permutation P
+        public Permutation P
         {
-            get
-            {
-                return Permutation.FromInversions(Pivots);
-            }
+            get { return _lazyP.Value; }
         }
 
         /// <summary>
@@ -108,12 +111,6 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization
         /// <returns>The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</returns>
         public virtual Matrix<T> Solve(Matrix<T> input)
         {
-            // Check for proper arguments.
-            if (input == null)
-            {
-                throw new ArgumentNullException("input");
-            }
-
             var x = input.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(input, x);
             return x;
@@ -133,12 +130,6 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization
         /// <returns>The left hand side <see cref="Vector{T}"/>, <b>x</b>.</returns>
         public virtual Vector<T> Solve(Vector<T> input)
         {
-            // Check for proper arguments.
-            if (input == null)
-            {
-                throw new ArgumentNullException("input");
-            }
-
             var x = input.CreateVector(input.Count);
             Solve(input, x);
             return x;
diff --git a/src/Numerics/LinearAlgebra/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Factorization/QR.cs
index 20803d56..d75fe570 100644
--- a/src/Numerics/LinearAlgebra/Factorization/QR.cs
+++ b/src/Numerics/LinearAlgebra/Factorization/QR.cs
@@ -60,6 +60,13 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization
     public abstract class QR<T> : ISolver<T>
     where T : struct, IEquatable<T>, IFormattable
     {
+        readonly Lazy<Matrix<T>> _lazyR;
+
+        protected QR()
+        {
+            _lazyR = new Lazy<Matrix<T>>(ComputeR);
+        }
+
         /// <summary>
         /// Gets or sets orthogonal Q matrix
         /// </summary>
@@ -80,7 +87,7 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization
         /// </summary>
         public Matrix<T> R
         {
-            get { return MatrixR.UpperTriangle(); }
+            get { return _lazyR.Value; }
         }
 
         /// <summary>
@@ -94,6 +101,11 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization
         /// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
         public abstract bool IsFullRank { get; }
 
+        private Matrix<T> ComputeR()
+        {
+            return MatrixR.UpperTriangle();
+        }
+
         /// <summary>
         /// Solves a system of linear equations, <b>AX = B</b>, with A QR factorized.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 78c27dd7..c6033e84 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -1042,7 +1042,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
         public override LU<float> LU()
         {
-            return new DenseLU(this);
+            return DenseLU.Create(this);
         }
 
         public override QR<float> QR(QRMethod method = QRMethod.Thin)
diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs
index 383053e9..46d256a9 100644
--- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs
+++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 {
@@ -41,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class DenseLU : LU
+    public sealed class DenseLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseLU"/> class. This object will compute the
@@ -50,7 +50,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// <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 DenseLU(DenseMatrix matrix)
+        public static DenseLU Create(DenseMatrix matrix)
         {
             if (matrix == null)
             {
@@ -63,12 +63,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             }
 
             // Create an array for the pivot indices.
-            Pivots = new int[matrix.RowCount];
+            var pivots = new int[matrix.RowCount];
 
             // Create a new matrix for the LU factors, then perform factorization (while overwriting).
-            var factors = (DenseMatrix)matrix.Clone();
-            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, Pivots);
-            Factors = factors;
+            var factors = (DenseMatrix) matrix.Clone();
+            Control.LinearAlgebraProvider.LUFactor(factors.Values, factors.RowCount, pivots);
+
+            return new DenseLU(factors, pivots);
+        }
+
+        DenseLU(Matrix<float> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -118,10 +124,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             }
 
             // Copy the contents of input to result.
-            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat);
+            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
+            var dfactors = (DenseMatrix) Factors;
             Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
@@ -167,10 +173,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             }
 
             // Copy the contents of input to result.
-            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat);
+            Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat);
 
             // LU solve by overwriting result.
-            var dfactors = (DenseMatrix)Factors;
+            var dfactors = (DenseMatrix) Factors;
             Control.LinearAlgebraProvider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
         }
 
@@ -180,7 +186,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// <returns>The inverse of this matrix.</returns>
         public override Matrix<float> Inverse()
         {
-            var result = (DenseMatrix)Factors.Clone();
+            var result = (DenseMatrix) Factors.Clone();
             Control.LinearAlgebraProvider.LUInverseFactored(result.Values, result.RowCount, Pivots);
             return result;
         }
diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/LU.cs
index 792417d8..dab327de 100644
--- a/src/Numerics/LinearAlgebra/Single/Factorization/LU.cs
+++ b/src/Numerics/LinearAlgebra/Single/Factorization/LU.cs
@@ -3,7 +3,9 @@
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
-// Copyright (c) 2009-2010 Math.NET
+//
+// Copyright (c) 2009-2013 Math.NET
+//
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
 // files (the "Software"), to deal in the Software without
@@ -12,8 +14,10 @@
 // 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
@@ -40,6 +44,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
     /// </remarks>
     public abstract class LU : LU<float>
     {
+        protected LU(Matrix<float> factors, int[] pivots)
+            : base(factors, pivots)
+        {
+        }
+
         /// <summary>
         /// Gets the determinant of the matrix for which the LU factorization was computed.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs
index ac250200..2689491e 100644
--- a/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs
+++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserLU.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2013 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -28,8 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-using MathNet.Numerics.Properties;
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 {
@@ -41,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
     /// <remarks>
     /// The computation of the LU factorization is done at construction time.
     /// </remarks>
-    public class UserLU : LU
+    public sealed class UserLU : LU
     {
         /// <summary>
         /// Initializes a new instance of the <see cref="UserLU"/> class. This object will compute the
@@ -50,7 +50,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// <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<float> matrix)
+        public static UserLU Create(Matrix<float> matrix)
         {
             if (matrix == null)
             {
@@ -64,13 +64,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 
             // Create an array for the pivot indices.
             var order = matrix.RowCount;
-            Factors = matrix.Clone();
-            Pivots = new int[order];
-            
+            var factors = matrix.Clone();
+            var pivots = new int[order];
+
             // Initialize the pivot matrix to the identity permutation.
             for (var i = 0; i < order; i++)
             {
-                Pivots[i] = i;
+                pivots[i] = i;
             }
 
             var vectorLUcolj = new float[order];
@@ -79,7 +79,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 // Make a copy of the j-th column to localize references.
                 for (var i = 0; i < order; i++)
                 {
-                    vectorLUcolj[i] = Factors.At(i, j);
+                    vectorLUcolj[i] = factors.At(i, j);
                 }
 
                 // Apply previous transformations.
@@ -89,11 +89,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                     var s = 0.0f;
                     for (var k = 0; k < kmax; k++)
                     {
-                        s += Factors.At(i, k) * vectorLUcolj[k];
+                        s += factors.At(i, k)*vectorLUcolj[k];
                     }
 
                     vectorLUcolj[i] -= s;
-                    Factors.At(i, j, vectorLUcolj[i]);
+                    factors.At(i, j, vectorLUcolj[i]);
                 }
 
                 // Find pivot and exchange if necessary.
@@ -110,23 +110,30 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 {
                     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);
+                        var temp = factors.At(p, k);
+                        factors.At(p, k, factors.At(j, k));
+                        factors.At(j, k, temp);
                     }
 
-                    Pivots[j] = p;
+                    pivots[j] = p;
                 }
 
                 // Compute multipliers.
-                if (j < order & Factors.At(j, j) != 0.0)
+                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)));
+                        factors.At(i, j, (factors.At(i, j)/factors.At(j, j)));
                     }
                 }
             }
+
+            return new UserLU(factors, pivots);
+        }
+
+        UserLU(Matrix<float> factors, int[] pivots)
+            : base(factors, pivots)
+        {
         }
 
         /// <summary>
@@ -182,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             }
 
             var order = Factors.RowCount;
-            
+
             // Solve L*Y = P*B
             for (var k = 0; k < order; k++)
             {
@@ -190,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 {
                     for (var j = 0; j < result.ColumnCount; j++)
                     {
-                        var temp = result.At(k, j) * Factors.At(i, k);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -201,14 +208,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             {
                 for (var j = 0; j < result.ColumnCount; j++)
                 {
-                    result.At(k, j, (result.At(k, j) / Factors.At(k, k)));
+                    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);
+                        var temp = result.At(k, j)*Factors.At(i, k);
                         result.At(i, j, result.At(i, j) - temp);
                     }
                 }
@@ -258,7 +265,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 result[p] = result[i];
                 result[i] = temp;
             }
-            
+
             var order = Factors.RowCount;
 
             // Solve L*Y = P*B
@@ -266,7 +273,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             {
                 for (var i = k + 1; i < order; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
 
@@ -276,7 +283,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 result[k] /= Factors.At(k, k);
                 for (var i = 0; i < k; i++)
                 {
-                    result[i] -= result[k] * Factors.At(i, k);
+                    result[i] -= result[k]*Factors.At(i, k);
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs
index c0d2ab08..5486e14f 100644
--- a/src/Numerics/LinearAlgebra/Single/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs
@@ -466,7 +466,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
         public override LU<float> LU()
         {
-            return new UserLU(this);
+            return UserLU.Create(this);
         }
 
         public override QR<float> QR(QRMethod method = QRMethod.Thin)