diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index 44e0f932..c1a47d7e 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -149,6 +149,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DenseMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DenseMatrix(Matrix<Complex> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index 7c211073..dd9d7885 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -149,6 +149,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DiagonalMatrix(Matrix<Complex> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
index 0ef336ae..48dd3926 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
@@ -184,32 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public SparseMatrix(Matrix<Complex> matrix)
             : this(matrix.RowCount, matrix.ColumnCount)
         {
-            var sparseMatrix = matrix as SparseMatrix;
-
-            var rows = matrix.RowCount;
-            var columns = matrix.ColumnCount;
-
-            if (sparseMatrix == null)
-            {
-                for (var i = 0; i < rows; i++)
-                {
-                    for (var j = 0; j < columns; j++)
-                    {
-                        _storage.At(i, j, matrix.At(i, j));
-                    }
-                }
-            }
-            else
-            {
-                var matrixStorage = sparseMatrix.Raw;
-                var valueCount = _storage.ValueCount = matrixStorage.ValueCount;
-                _storage.ColumnIndices = new int[valueCount];
-                _storage.Values = new Complex[valueCount];
-
-                Array.Copy(matrixStorage.Values, _storage.Values, valueCount);
-                Array.Copy(matrixStorage.ColumnIndices, _storage.ColumnIndices, valueCount);
-                Array.Copy(matrixStorage.RowPointers, _storage.RowPointers, rows);
-            }
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
         }
 
         /// <summary>
@@ -645,10 +620,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseMatrix(ColumnCount, RowCount);
-            var retStorage = ret.Raw;
-            retStorage.ColumnIndices = new int[valueCount];
-            retStorage.Values = new Complex[valueCount];
+            var ret = new SparseCompressedRowMatrixStorage<Complex>(ColumnCount, RowCount, Complex.Zero)
+                {
+                    ColumnIndices = new int[valueCount],
+                    Values = new Complex[valueCount]
+                };
 
             // Do an 'inverse' CopyTo iterate over the rows
             for (var i = 0; i < rowPointers.Length; i++)
@@ -666,20 +642,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
                 for (var j = startIndex; j < endIndex; j++)
                 {
-                    retStorage.At(columnIndices[j], i, values[j]);
+                    ret.At(columnIndices[j], i, values[j]);
                 }
             }
 
-            return ret;
+            return new SparseMatrix(ret);
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override Complex FrobeniusNorm()
         {
-            var transpose = (SparseMatrix)Transpose();
-            var aat = (this * transpose).Raw;
-
+            var aat = (SparseCompressedRowMatrixStorage<Complex>) (this*Transpose()).Storage;
             var norm = 0d;
 
             for (var i = 0; i < aat.RowPointers.Length; i++)
@@ -875,21 +849,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var m = new SparseMatrix(order);
-            var mStorage = m.Raw;
-
-            mStorage.ValueCount = order;
-            mStorage.Values = new Complex[order];
-            mStorage.ColumnIndices = new int[order];
+            var m = new SparseCompressedRowMatrixStorage<Complex>(order, order, Complex.Zero)
+                {
+                    ValueCount = order,
+                    Values = new Complex[order],
+                    ColumnIndices = new int[order]
+                };
 
             for (var i = 0; i < order; i++)
             {
-                mStorage.Values[i] = 1d;
-                mStorage.ColumnIndices[i] = i;
-                mStorage.RowPointers[i] = i;
+                m.Values[i] = 1d;
+                m.ColumnIndices[i] = i;
+                m.RowPointers[i] = i;
             }
 
-            return m;
+            return new SparseMatrix(m);
         }
         #endregion
 
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 66309118..56ed90ac 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -149,6 +149,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DenseMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DenseMatrix(Matrix<Complex32> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index 8160fed5..fac1a33b 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
@@ -149,6 +149,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DiagonalMatrix(Matrix<Complex32> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
index 74244ed8..be05bce5 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
@@ -184,32 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public SparseMatrix(Matrix<Complex32> matrix)
             : this(matrix.RowCount, matrix.ColumnCount)
         {
-            var sparseMatrix = matrix as SparseMatrix;
-
-            var rows = matrix.RowCount;
-            var columns = matrix.ColumnCount;
-
-            if (sparseMatrix == null)
-            {
-                for (var i = 0; i < rows; i++)
-                {
-                    for (var j = 0; j < columns; j++)
-                    {
-                        _storage.At(i, j, matrix.At(i, j));
-                    }
-                }
-            }
-            else
-            {
-                var matrixStorage = sparseMatrix.Raw;
-                var valueCount = _storage.ValueCount = matrixStorage.ValueCount;
-                _storage.ColumnIndices = new int[valueCount];
-                _storage.Values = new Complex32[valueCount];
-
-                Array.Copy(matrixStorage.Values, _storage.Values, valueCount);
-                Array.Copy(matrixStorage.ColumnIndices, _storage.ColumnIndices, valueCount);
-                Array.Copy(matrixStorage.RowPointers, _storage.RowPointers, rows);
-            }
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
         }
 
         /// <summary>
@@ -645,10 +620,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseMatrix(ColumnCount, RowCount);
-            var retStorage = ret.Raw;
-            retStorage.ColumnIndices = new int[valueCount];
-            retStorage.Values = new Complex32[valueCount];
+            var ret = new SparseCompressedRowMatrixStorage<Complex32>(ColumnCount, RowCount, Complex32.Zero)
+                {
+                    ColumnIndices = new int[valueCount],
+                    Values = new Complex32[valueCount]
+                };
 
             // Do an 'inverse' CopyTo iterate over the rows
             for (var i = 0; i < rowPointers.Length; i++)
@@ -666,20 +642,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
                 for (var j = startIndex; j < endIndex; j++)
                 {
-                    retStorage.At(columnIndices[j], i, values[j]);
+                    ret.At(columnIndices[j], i, values[j]);
                 }
             }
 
-            return ret;
+            return new SparseMatrix(ret);
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override Complex32 FrobeniusNorm()
         {
-            var transpose = (SparseMatrix)Transpose();
-            var aat = (this * transpose).Raw;
-
+            var aat = (SparseCompressedRowMatrixStorage<Complex32>) (this*Transpose()).Storage;
             var norm = 0f;
 
             for (var i = 0; i < aat.RowPointers.Length; i++)
@@ -875,21 +849,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var m = new SparseMatrix(order);
-            var mStorage = m.Raw;
-
-            mStorage.ValueCount = order;
-            mStorage.Values = new Complex32[order];
-            mStorage.ColumnIndices = new int[order];
+            var m = new SparseCompressedRowMatrixStorage<Complex32>(order, order, Complex32.Zero)
+                {
+                    ValueCount = order,
+                    Values = new Complex32[order],
+                    ColumnIndices = new int[order]
+                };
 
             for (var i = 0; i < order; i++)
             {
-                mStorage.Values[i] = 1f;
-                mStorage.ColumnIndices[i] = i;
-                mStorage.RowPointers[i] = i;
+                m.Values[i] = 1f;
+                m.ColumnIndices[i] = i;
+                m.RowPointers[i] = i;
             }
 
-            return m;
+            return new SparseMatrix(m);
         }
         #endregion
 
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index 71fff598..a965b53a 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -149,6 +149,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DenseMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DenseMatrix(Matrix<double> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index ad914a0d..f43217f6 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -148,6 +148,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DiagonalMatrix(Matrix<double> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index f5a5bd82..1e702faa 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
@@ -183,32 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public SparseMatrix(Matrix<double> matrix)
             : this(matrix.RowCount, matrix.ColumnCount)
         {
-            var sparseMatrix = matrix as SparseMatrix;
-
-            var rows = matrix.RowCount;
-            var columns = matrix.ColumnCount;
-
-            if (sparseMatrix == null)
-            {
-                for (var i = 0; i < rows; i++)
-                {
-                    for (var j = 0; j < columns; j++)
-                    {
-                        _storage.At(i, j, matrix.At(i, j));
-                    }
-                }  
-            }
-            else
-            {
-                var matrixStorage = sparseMatrix.Raw;
-                var valueCount = _storage.ValueCount = matrixStorage.ValueCount;
-                _storage.ColumnIndices = new int[valueCount];
-                _storage.Values = new double[valueCount];
-
-                Buffer.BlockCopy(matrixStorage.Values, 0, _storage.Values, 0, valueCount * Constants.SizeOfDouble);
-                Buffer.BlockCopy(matrixStorage.ColumnIndices, 0, _storage.ColumnIndices, 0, valueCount * Constants.SizeOfInt);
-                Buffer.BlockCopy(matrixStorage.RowPointers, 0, _storage.RowPointers, 0, rows * Constants.SizeOfInt);
-            }
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
         }
 
         /// <summary>
@@ -644,10 +619,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseMatrix(ColumnCount, RowCount);
-            var retStorage = ret.Raw;
-            retStorage.ColumnIndices = new int[valueCount];
-            retStorage.Values = new double[valueCount];
+            var ret = new SparseCompressedRowMatrixStorage<double>(ColumnCount, RowCount, 0d)
+                {
+                    ColumnIndices = new int[valueCount],
+                    Values = new double[valueCount]
+                };
 
             // Do an 'inverse' CopyTo iterate over the rows
             for (var i = 0; i < rowPointers.Length; i++)
@@ -665,20 +641,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
                 for (var j = startIndex; j < endIndex; j++)
                 {
-                    retStorage.At(columnIndices[j], i, values[j]);
+                    ret.At(columnIndices[j], i, values[j]);
                 }
             }
 
-            return ret;
+            return new SparseMatrix(ret);
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override double FrobeniusNorm()
         {
-            var transpose = (SparseMatrix)Transpose();
-            var aat = (this * transpose).Raw;
-
+            var aat = (SparseCompressedRowMatrixStorage<double>) (this*Transpose()).Storage;
             var norm = 0d;
 
             for (var i = 0; i < aat.RowPointers.Length; i++)
@@ -873,21 +847,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var m = new SparseMatrix(order);
-            var mStorage = m.Raw;
-
-            mStorage.ValueCount = order;
-            mStorage.Values = new double[order];
-            mStorage.ColumnIndices = new int[order];
+            var m = new SparseCompressedRowMatrixStorage<double>(order, order, 0d)
+                {
+                    ValueCount = order,
+                    Values = new double[order],
+                    ColumnIndices = new int[order]
+                };
 
             for (var i = 0; i < order; i++)
             {
-                mStorage.Values[i] = 1d;
-                mStorage.ColumnIndices[i] = i;
-                mStorage.RowPointers[i] = i;
+                m.Values[i] = 1d;
+                m.ColumnIndices[i] = i;
+                m.RowPointers[i] = i;
             }
 
-            return m;
+            return new SparseMatrix(m);
         }
         #endregion
 
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index cd458836..ea428a0c 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -149,6 +149,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DenseMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DenseMatrix(Matrix<float> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index 3dc10eaa..41d84eda 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -148,6 +148,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class, copying
+        /// the values from the given matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to copy.</param>
+        public DiagonalMatrix(Matrix<float> matrix)
+            : this(matrix.RowCount, matrix.ColumnCount)
+        {
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
+        }
+
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index 4c89cc1e..60719210 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
@@ -183,32 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public SparseMatrix(Matrix<float> matrix)
             : this(matrix.RowCount, matrix.ColumnCount)
         {
-            var sparseMatrix = matrix as SparseMatrix;
-
-            var rows = matrix.RowCount;
-            var columns = matrix.ColumnCount;
-
-            if (sparseMatrix == null)
-            {
-                for (var i = 0; i < rows; i++)
-                {
-                    for (var j = 0; j < columns; j++)
-                    {
-                        _storage.At(i, j, matrix.At(i, j));
-                    }
-                }
-            }
-            else
-            {
-                var matrixStorage = sparseMatrix.Raw;
-                var valueCount = _storage.ValueCount = matrixStorage.ValueCount;
-                _storage.ColumnIndices = new int[valueCount];
-                _storage.Values = new float[valueCount];
-
-                Buffer.BlockCopy(matrixStorage.Values, 0, _storage.Values, 0, valueCount * Constants.SizeOfFloat);
-                Buffer.BlockCopy(matrixStorage.ColumnIndices, 0, _storage.ColumnIndices, 0, valueCount * Constants.SizeOfInt);
-                Buffer.BlockCopy(matrixStorage.RowPointers, 0, _storage.RowPointers, 0, rows * Constants.SizeOfInt);
-            }
+            matrix.Storage.CopyTo(Storage, skipClearing: true);
         }
 
         /// <summary>
@@ -644,10 +619,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseMatrix(ColumnCount, RowCount);
-            var retStorage = ret.Raw;
-            retStorage.ColumnIndices = new int[valueCount];
-            retStorage.Values = new float[valueCount];
+            var ret = new SparseCompressedRowMatrixStorage<float>(ColumnCount, RowCount, 0f)
+                {
+                    ColumnIndices = new int[valueCount],
+                    Values = new float[valueCount]
+                };
 
             // Do an 'inverse' CopyTo iterate over the rows
             for (var i = 0; i < rowPointers.Length; i++)
@@ -665,20 +641,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
                 for (var j = startIndex; j < endIndex; j++)
                 {
-                    retStorage.At(columnIndices[j], i, values[j]);
+                    ret.At(columnIndices[j], i, values[j]);
                 }
             }
 
-            return ret;
+            return new SparseMatrix(ret);
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override float FrobeniusNorm()
         {
-            var transpose = (SparseMatrix)Transpose();
-            var aat = (this * transpose).Raw;
-
+            var aat = (SparseCompressedRowMatrixStorage<float>) (this*Transpose()).Storage;
             var norm = 0f;
 
             for (var i = 0; i < aat.RowPointers.Length; i++)
@@ -873,12 +847,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var m = new SparseMatrix(order);
-            var mStorage = m.Raw;
-
-            mStorage.ValueCount = order;
-            mStorage.Values = new float[order];
-            mStorage.ColumnIndices = new int[order];
+            var mStorage = new SparseCompressedRowMatrixStorage<float>(order, order, 0f)
+                {
+                    ValueCount = order,
+                    Values = new float[order],
+                    ColumnIndices = new int[order]
+                };
 
             for (var i = 0; i < order; i++)
             {
@@ -887,7 +861,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 mStorage.RowPointers[i] = i;
             }
 
-            return m;
+            return new SparseMatrix(mStorage);
         }
         #endregion