diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index c1a47d7e..57f10f58 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -198,59 +198,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new DenseVector(size);
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<Complex> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<Complex>(rowCount, columnCount);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount);
-            return new DenseMatrix(storage.RowCount, storage.ColumnCount, storage.Data);
-        }
-
-        /// <summary>
-        /// Copies the values of a given matrix into a region in this matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying to.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying to.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <param name="subMatrix">The sub-matrix to copy from.</param>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception> 
-        /// <exception cref="ArgumentNullException">If <paramref name="subMatrix"/> is <see langword="null" /></exception>
-        /// <item>the size of <paramref name="subMatrix"/> is not at least <paramref name="rowCount"/> x <paramref name="columnCount"/>.</item>
-        /// <exception cref="ArgumentException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<Complex> subMatrix)
-        {
-            var denseSubMatrix = subMatrix as DenseMatrix;
-            if (denseSubMatrix != null)
-            {
-                denseSubMatrix._storage.CopySubMatrixTo(_storage, 0, rowIndex, rowCount, 0, columnIndex, columnCount);
-                return;
-            }
-
-            base.SetSubMatrix(rowIndex, rowCount, columnIndex, columnCount, subMatrix);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>        
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index d463e5d8..214a5a87 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -1027,11 +1027,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// is not positive.</exception>
         public override Matrix<Complex> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
+            var target = rowIndex == columnIndex
+                ? (Matrix<Complex>)new DiagonalMatrix(rowCount, columnCount)
+                : new SparseMatrix(rowCount, columnCount);
 
-            var storage = new SparseCompressedRowMatrixStorage<Complex>(rowCount, columnCount, Complex.Zero);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-            return new SparseMatrix(storage);
+            Storage.CopySubMatrixTo(target.Storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, skipClearing: true);
+            return target;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
index 4dd867fb..6e7275a2 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
@@ -350,83 +350,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<Complex> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            if (rowIndex >= RowCount || rowIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("rowCount", Resources.ArgumentMustBePositive);
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("columnCount", Resources.ArgumentMustBePositive);
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (colMax > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            var result = (SparseMatrix)CreateMatrix(rowCount, columnCount);
-
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (int i = rowIndex, row = 0; i < rowMax; i++, row++)
-            {
-                var startIndex = rowPointers[i];
-                var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
-
-                for (int j = startIndex; j < endIndex; j++)
-                {
-                    // check if the column index is in the range
-                    if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount))
-                    {
-                        var column = columnIndices[j] - columnIndex;
-                        result._storage.At(row, column, values[j]);
-                    }
-                }
-            }
-
-            return result;
-        }
-
         /// <summary>
         /// Returns a new matrix containing the lower triangle of this matrix. The new matrix
         /// does not contain the diagonal elements of this matrix.
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 56ed90ac..8fe6314f 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -198,59 +198,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new DenseVector(size);
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<Complex32> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<Complex32>(rowCount, columnCount);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount);
-            return new DenseMatrix(storage.RowCount, storage.ColumnCount, storage.Data);
-        }
-
-        /// <summary>
-        /// Copies the values of a given matrix into a region in this matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying to.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying to.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <param name="subMatrix">The sub-matrix to copy from.</param>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception> 
-        /// <exception cref="ArgumentNullException">If <paramref name="subMatrix"/> is <see langword="null" /></exception>
-        /// <item>the size of <paramref name="subMatrix"/> is not at least <paramref name="rowCount"/> x <paramref name="columnCount"/>.</item>
-        /// <exception cref="ArgumentException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<Complex32> subMatrix)
-        {
-            var denseSubMatrix = subMatrix as DenseMatrix;
-            if (denseSubMatrix != null)
-            {
-                denseSubMatrix._storage.CopySubMatrixTo(_storage, 0, rowIndex, rowCount, 0, columnIndex, columnCount);
-                return;
-            }
-
-            base.SetSubMatrix(rowIndex, rowCount, columnIndex, columnCount, subMatrix);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>        
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index 34fa4184..400d3f50 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
@@ -1027,11 +1027,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// is not positive.</exception>
         public override Matrix<Complex32> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
+            var target = rowIndex == columnIndex
+                ? (Matrix<Complex32>)new DiagonalMatrix(rowCount, columnCount)
+                : new SparseMatrix(rowCount, columnCount);
 
-            var storage = new SparseCompressedRowMatrixStorage<Complex32>(rowCount, columnCount, Complex32.Zero);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-            return new SparseMatrix(storage);
+            Storage.CopySubMatrixTo(target.Storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, skipClearing: true);
+            return target;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
index 224b7b2c..8164cb4b 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
@@ -350,83 +350,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<Complex32> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            if (rowIndex >= RowCount || rowIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("rowCount", Resources.ArgumentMustBePositive);
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("columnCount", Resources.ArgumentMustBePositive);
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (colMax > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            var result = (SparseMatrix)CreateMatrix(rowCount, columnCount);
-
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (int i = rowIndex, row = 0; i < rowMax; i++, row++)
-            {
-                var startIndex = rowPointers[i];
-                var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
-
-                for (int j = startIndex; j < endIndex; j++)
-                {
-                    // check if the column index is in the range
-                    if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount))
-                    {
-                        var column = columnIndices[j] - columnIndex;
-                        result._storage.At(row, column, values[j]);
-                    }
-                }
-            }
-
-            return result;
-        }
-
         /// <summary>
         /// Returns a new matrix containing the lower triangle of this matrix. The new matrix
         /// does not contain the diagonal elements of this matrix.
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index a965b53a..16756af5 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -198,59 +198,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new DenseVector(size);
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<double> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<double>(rowCount, columnCount);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount);
-            return new DenseMatrix(storage.RowCount, storage.ColumnCount, storage.Data);
-        }
-
-        /// <summary>
-        /// Copies the values of a given matrix into a region in this matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying to.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying to.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <param name="subMatrix">The sub-matrix to copy from.</param>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception> 
-        /// <exception cref="ArgumentNullException">If <paramref name="subMatrix"/> is <see langword="null" /></exception>
-        /// <item>the size of <paramref name="subMatrix"/> is not at least <paramref name="rowCount"/> x <paramref name="columnCount"/>.</item>
-        /// <exception cref="ArgumentException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<double> subMatrix)
-        {
-            var denseSubMatrix = subMatrix as DenseMatrix;
-            if (denseSubMatrix != null)
-            {
-                denseSubMatrix._storage.CopySubMatrixTo(_storage, 0, rowIndex, rowCount, 0, columnIndex, columnCount);
-                return;
-            }
-
-            base.SetSubMatrix(rowIndex, rowCount, columnIndex, columnCount, subMatrix);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>        
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 5832e075..4ef3b400 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -1021,11 +1021,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// is not positive.</exception>
         public override Matrix<double> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
+            var target = rowIndex == columnIndex
+                ? (Matrix<double>)new DiagonalMatrix(rowCount, columnCount)
+                : new SparseMatrix(rowCount, columnCount);
 
-            var storage = new SparseCompressedRowMatrixStorage<double>(rowCount, columnCount, 0d);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-            return new SparseMatrix(storage);
+            Storage.CopySubMatrixTo(target.Storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, skipClearing: true);
+            return target;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index dc7ead37..8b9cf66a 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
@@ -349,83 +349,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<double> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            if (rowIndex >= RowCount || rowIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("rowCount", Resources.ArgumentMustBePositive);
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("columnCount", Resources.ArgumentMustBePositive);
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (colMax > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            var result = (SparseMatrix)CreateMatrix(rowCount, columnCount);
-
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (int i = rowIndex, row = 0; i < rowMax; i++, row++)
-            {
-                var startIndex = rowPointers[i];
-                var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
-
-                for (int j = startIndex; j < endIndex; j++)
-                {
-                    // check if the column index is in the range
-                    if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount))
-                    {
-                        var column = columnIndices[j] - columnIndex;
-                        result._storage.At(row, column, values[j]);
-                    }
-                }
-            }
-
-            return result;
-        }
-
         /// <summary>
         /// Returns a new matrix containing the lower triangle of this matrix. The new matrix
         /// does not contain the diagonal elements of this matrix.
diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.cs
index eeb4edb4..0c230c1a 100644
--- a/src/Numerics/LinearAlgebra/Generic/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Generic/Matrix.cs
@@ -221,7 +221,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         public virtual Matrix<T> Clone()
         {
             var result = CreateMatrix(RowCount, ColumnCount);
-            Storage.CopyTo(result.Storage);
+            Storage.CopyTo(result.Storage, skipClearing: true);
             return result;
         }
 
@@ -615,50 +615,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// is not positive.</exception>
         public virtual Matrix<T> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            if (rowIndex >= RowCount || rowIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("rowCount", Resources.ArgumentMustBePositive);
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("columnCount", Resources.ArgumentMustBePositive);
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (colMax > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            var result = CreateMatrix(rowCount, columnCount);
-
-            for (var j = columnIndex; j < colMax; j++)
-            {
-                for (int i = rowIndex, ii = 0; i < rowMax; i++, ii++)
-                {
-                    result.At(ii, j - columnIndex, At(i, j));
-                }
-            }
-
-            return result;
+            var target = CreateMatrix(rowCount, columnCount);
+            Storage.CopySubMatrixTo(target.Storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount);
+            return target;
         }
 
         /// <summary>
@@ -1113,61 +1072,12 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// is not positive.</exception>
         public virtual void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<T> subMatrix)
         {
-            if (rowIndex >= RowCount || rowIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("rowCount", Resources.ArgumentMustBePositive);
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("columnCount", Resources.ArgumentMustBePositive);
-            }
-
             if (subMatrix == null)
             {
                 throw new ArgumentNullException("subMatrix");
             }
 
-            if (columnCount > subMatrix.ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount", @"columnLength can be at most the number of columns in subMatrix.");
-            }
-
-            if (rowCount > subMatrix.RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount", @"rowLength can be at most the number of rows in subMatrix.");
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (colMax > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            for (var j = columnIndex; j < colMax; j++)
-            {
-                for (int i = rowIndex, ii = 0; i < rowMax; i++, ii++)
-                {
-                    At(i, j, subMatrix.At(ii, j - columnIndex));
-                }
-            }
+            subMatrix.Storage.CopySubMatrixTo(Storage, 0, rowIndex, rowCount, 0, columnIndex, columnCount);
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index ea428a0c..2e4c9463 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -198,59 +198,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new DenseVector(size);
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<float> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<float>(rowCount, columnCount);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount);
-            return new DenseMatrix(storage.RowCount, storage.ColumnCount, storage.Data);
-        }
-
-        /// <summary>
-        /// Copies the values of a given matrix into a region in this matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying to.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying to.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <param name="subMatrix">The sub-matrix to copy from.</param>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception> 
-        /// <exception cref="ArgumentNullException">If <paramref name="subMatrix"/> is <see langword="null" /></exception>
-        /// <item>the size of <paramref name="subMatrix"/> is not at least <paramref name="rowCount"/> x <paramref name="columnCount"/>.</item>
-        /// <exception cref="ArgumentException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<float> subMatrix)
-        {
-            var denseSubMatrix = subMatrix as DenseMatrix;
-            if (denseSubMatrix != null)
-            {
-                denseSubMatrix._storage.CopySubMatrixTo(_storage, 0, rowIndex, rowCount, 0, columnIndex, columnCount);
-                return;
-            }
-
-            base.SetSubMatrix(rowIndex, rowCount, columnIndex, columnCount, subMatrix);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>        
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index f43e28b5..aa15c269 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -1021,11 +1021,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// is not positive.</exception>
         public override Matrix<float> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
+            var target = rowIndex == columnIndex
+                ? (Matrix<float>)new DiagonalMatrix(rowCount, columnCount)
+                : new SparseMatrix(rowCount, columnCount);
 
-            var storage = new SparseCompressedRowMatrixStorage<float>(rowCount, columnCount, 0f);
-            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-            return new SparseMatrix(storage);
+            Storage.CopySubMatrixTo(target.Storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, skipClearing: true);
+            return target;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index 7aa3bcfc..4f1d6b5a 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
@@ -349,83 +349,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
-        /// <summary>
-        /// Creates a matrix that contains the values from the requested sub-matrix.
-        /// </summary>
-        /// <param name="rowIndex">The row to start copying from.</param>
-        /// <param name="rowCount">The number of rows to copy. Must be positive.</param>
-        /// <param name="columnIndex">The column to start copying from.</param>
-        /// <param name="columnCount">The number of columns to copy. Must be positive.</param>
-        /// <returns>The requested sub-matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
-        /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
-        /// of columns.</item>
-        /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
-        /// is not positive.</exception>
-        public override Matrix<float> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
-        {
-            if (rowIndex >= RowCount || rowIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("rowCount", Resources.ArgumentMustBePositive);
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentOutOfRangeException("columnCount", Resources.ArgumentMustBePositive);
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (colMax > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            var result = (SparseMatrix)CreateMatrix(rowCount, columnCount);
-
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (int i = rowIndex, row = 0; i < rowMax; i++, row++)
-            {
-                var startIndex = rowPointers[i];
-                var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
-
-                for (int j = startIndex; j < endIndex; j++)
-                {
-                    // check if the column index is in the range
-                    if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount))
-                    {
-                        var column = columnIndices[j] - columnIndex;
-                        result._storage.At(row, column, values[j]);
-                    }
-                }
-            }
-
-            return result;
-        }
-
         /// <summary>
         /// Returns a new matrix containing the lower triangle of this matrix. The new matrix
         /// does not contain the diagonal elements of this matrix.
diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
index 45f29a4f..e21c503b 100644
--- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
@@ -53,6 +53,20 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Clear(Data, 0, Data.Length);
         }
 
+        public override void Clear(int rowIndex, int rowCount, int columnIndex, int columnCount)
+        {
+            if (rowIndex == 0 && columnIndex == 0 && rowCount == RowCount && columnCount == ColumnCount)
+            {
+                Clear();
+                return;
+            }
+
+            for (int j = columnIndex; j < columnIndex + columnCount; j++)
+            {
+                Array.Clear(Data, j*RowCount + rowIndex, rowCount);
+            }
+        }
+
         /// <remarks>Parameters assumed to be validated already.</remarks>
         public override void CopyTo(MatrixStorage<T> target, bool skipClearing = false)
         {
@@ -74,8 +88,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-
-        public void CopyTo(DenseColumnMajorMatrixStorage<T> target)
+        void CopyTo(DenseColumnMajorMatrixStorage<T> target)
         {
             if (ReferenceEquals(this, target))
             {
@@ -97,7 +110,24 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Copy(Data, 0, target.Data, 0, Data.Length);
         }
 
-        public void CopySubMatrixTo(DenseColumnMajorMatrixStorage<T> target,
+        public override void CopySubMatrixTo(MatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount,
+            bool skipClearing = false)
+        {
+            var denseTarget = target as DenseColumnMajorMatrixStorage<T>;
+            if (denseTarget != null)
+            {
+                CopySubMatrixTo(denseTarget, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount);
+                return;
+            }
+
+            // FALL BACK
+
+            base.CopySubMatrixTo(target, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing);
+        }
+
+        void CopySubMatrixTo(DenseColumnMajorMatrixStorage<T> target,
             int sourceRowIndex, int targetRowIndex, int rowCount,
             int sourceColumnIndex, int targetColumnIndex, int columnCount)
         {
diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
index 55b67bcd..cdd1f753 100644
--- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
@@ -74,6 +74,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Clear(Data, 0, Data.Length);
         }
 
+        public override void Clear(int rowIndex, int rowCount, int columnIndex, int columnCount)
+        {
+            var beginInclusive = Math.Max(rowIndex, columnIndex);
+            var endExclusive = Math.Min(rowIndex + rowCount, columnIndex + columnCount);
+            if (endExclusive > beginInclusive)
+            {
+                Array.Clear(Data, beginInclusive, endExclusive - beginInclusive);
+            }
+        }
+
         /// <summary>
         /// Indicates whether the current object is equal to another object of the same type.
         /// </summary>
@@ -169,7 +179,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        public void CopyTo(DiagonalMatrixStorage<T> target)
+        void CopyTo(DiagonalMatrixStorage<T> target)
         {
             if (ReferenceEquals(this, target))
             {
@@ -191,7 +201,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Copy(Data, 0, target.Data, 0, Data.Length);
         }
 
-        public void CopyTo(SparseCompressedRowMatrixStorage<T> target, bool skipClearing = false)
+        void CopyTo(SparseCompressedRowMatrixStorage<T> target, bool skipClearing = false)
         {
             if (target == null)
             {
@@ -215,7 +225,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        public void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing = false)
+        void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing = false)
         {
             if (target == null)
             {
@@ -239,7 +249,76 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        public void CopySubMatrixTo(DenseColumnMajorMatrixStorage<T> target,
+        public override void CopySubMatrixTo(MatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount,
+            bool skipClearing = false)
+        {
+            var denseTarget = target as DenseColumnMajorMatrixStorage<T>;
+            if (denseTarget != null)
+            {
+                CopySubMatrixTo(denseTarget, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing);
+                return;
+            }
+
+            var diagonalTarget = target as DiagonalMatrixStorage<T>;
+            if (diagonalTarget != null)
+            {
+                CopySubMatrixTo(diagonalTarget, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount);
+                return;
+            }
+
+            var sparseTarget = target as SparseCompressedRowMatrixStorage<T>;
+            if (sparseTarget != null)
+            {
+                CopySubMatrixTo(sparseTarget, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing);
+                return;
+            }
+
+            // FALL BACK
+
+            base.CopySubMatrixTo(target, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing);
+        }
+
+        void CopySubMatrixTo(DiagonalMatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (ReferenceEquals(this, target))
+            {
+                throw new NotSupportedException();
+            }
+
+            if (sourceRowIndex - sourceColumnIndex != targetRowIndex - targetColumnIndex)
+            {
+                if (Data.Any(x => !_zero.Equals(x)))
+                {
+                    throw new NotSupportedException();
+                }
+
+                target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
+                return;
+            }
+
+            ValidateSubMatrixRange(target,
+                sourceRowIndex, targetRowIndex, rowCount,
+                sourceColumnIndex, targetColumnIndex, columnCount);
+
+            var beginInclusive = Math.Max(sourceRowIndex, sourceColumnIndex);
+            var endExclusive = Math.Min(sourceRowIndex + rowCount, sourceColumnIndex + columnCount);
+            if (endExclusive > beginInclusive)
+            {
+                var beginTarget = Math.Max(targetRowIndex, targetColumnIndex);
+                Array.Copy(Data, beginInclusive, target.Data, beginTarget, endExclusive - beginInclusive);
+            }
+        }
+
+        void CopySubMatrixTo(DenseColumnMajorMatrixStorage<T> target,
             int sourceRowIndex, int targetRowIndex, int rowCount,
             int sourceColumnIndex, int targetColumnIndex, int columnCount,
             bool skipClearing = false)
@@ -255,7 +334,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             if (!skipClearing)
             {
-                target.Clear();
+                target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
             }
 
             if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex)
@@ -299,7 +378,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        public void CopySubMatrixTo(SparseCompressedRowMatrixStorage<T> target,
+        void CopySubMatrixTo(SparseCompressedRowMatrixStorage<T> target,
             int sourceRowIndex, int targetRowIndex, int rowCount,
             int sourceColumnIndex, int targetColumnIndex, int columnCount,
             bool skipClearing = false)
@@ -315,7 +394,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             if (!skipClearing)
             {
-                target.Clear();
+                target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
             }
 
             if (sourceRowIndex == sourceColumnIndex)
diff --git a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs
index b3ac5bab..a8dcbba1 100644
--- a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs
@@ -107,6 +107,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
+        public virtual void Clear(int rowIndex, int rowCount, int columnIndex, int columnCount)
+        {
+            for (var i = rowIndex; i < rowIndex + rowCount; i++)
+            {
+                for (var j = columnIndex; j < columnIndex + columnCount; j++)
+                {
+                    At(i, j, default(T));
+                }
+            }
+        }
+
         /// <summary>
         /// Indicates whether the current object is equal to another object of the same type.
         /// </summary>
@@ -194,5 +205,33 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 }
             }
         }
+
+        public virtual void CopySubMatrixTo(MatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount,
+            bool skipClearing = false)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (ReferenceEquals(this, target))
+            {
+                throw new NotSupportedException();
+            }
+
+            ValidateSubMatrixRange(target,
+                sourceRowIndex, targetRowIndex, rowCount,
+                sourceColumnIndex, targetColumnIndex, columnCount);
+
+            for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++)
+            {
+                for (int i = sourceRowIndex, ii = targetRowIndex; i < sourceRowIndex + rowCount; i++, ii++)
+                {
+                    target.At(ii, jj, At(i, j));
+                }
+            }
+        }
     }
 }
diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
index c7ac806b..ee0bbea2 100644
--- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
@@ -142,6 +142,33 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Clear(RowPointers, 0, RowPointers.Length);
         }
 
+        public override void Clear(int rowIndex, int rowCount, int columnIndex, int columnCount)
+        {
+            if (rowIndex == 0 && columnIndex == 0 && rowCount == RowCount && columnCount == ColumnCount)
+            {
+                Clear();
+                return;
+            }
+
+            for (int i = rowIndex, row = 0; i < rowIndex + rowCount; i++, row++)
+            {
+                var startIndex = RowPointers[i];
+                var endIndex = i < RowPointers.Length - 1 ? RowPointers[i + 1] : ValueCount;
+
+                for (int j = startIndex; j < endIndex; j++)
+                {
+                    // check if the column index is in the range
+                    if ((ColumnIndices[j] >= columnIndex) && (ColumnIndices[j] < columnIndex + columnCount))
+                    {
+                        var column = ColumnIndices[j] - columnIndex;
+
+                        // NOTE: potential for more efficient implementation
+                        At(row, column, _zero);
+                    }
+                }
+            }
+        }
+
         /// <summary>
         /// Delete value from internal storage
         /// </summary>
@@ -328,7 +355,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        public void CopyTo(SparseCompressedRowMatrixStorage<T> target)
+        void CopyTo(SparseCompressedRowMatrixStorage<T> target)
         {
             if (ReferenceEquals(this, target))
             {
@@ -358,7 +385,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        public void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing = false)
+        void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing = false)
         {
             if (target == null)
             {
@@ -389,5 +416,49 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 }
             }
         }
+
+        public override void CopySubMatrixTo(MatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount,
+            bool skipClearing = false)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (ReferenceEquals(this, target))
+            {
+                throw new NotSupportedException();
+            }
+
+            ValidateSubMatrixRange(target,
+                sourceRowIndex, targetRowIndex, rowCount,
+                sourceColumnIndex, targetColumnIndex, columnCount);
+
+            // NOTE: potential for more efficient implementation (specialized on target storage schema)
+            // (this would then essentially be the fallback implementation)
+
+            if (!skipClearing)
+            {
+                target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
+            }
+
+            for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++)
+            {
+                var startIndex = RowPointers[i];
+                var endIndex = i < RowPointers.Length - 1 ? RowPointers[i + 1] : ValueCount;
+
+                for (int j = startIndex; j < endIndex; j++)
+                {
+                    // check if the column index is in the range
+                    if ((ColumnIndices[j] >= sourceColumnIndex) && (ColumnIndices[j] < sourceColumnIndex + columnCount))
+                    {
+                        var column = ColumnIndices[j] - sourceColumnIndex;
+                        target.At(targetRowIndex + row, targetColumnIndex + column, Values[j]);
+                    }
+                }
+            }
+        }
     }
 }