diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index 57f10f58..e1afd250 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -172,16 +172,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <summary>
         /// Creates a <c>DenseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DenseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new DenseMatrix(numberOfRows, numberOfColumns);
         }
@@ -190,10 +187,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// Creates a <see cref="Vector{T}"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<Complex> CreateVector(int size)
+        public override Vector<Complex> CreateVector(int size, bool fullyMutable = false)
         {
             return new DenseVector(size);
         }
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index 214a5a87..805720e7 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -163,18 +163,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DiagonalMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
-            return new DiagonalMatrix(numberOfRows, numberOfColumns);
+            return fullyMutable
+                ? (Matrix<Complex>) new SparseMatrix(numberOfRows, numberOfColumns)
+                : new DiagonalMatrix(numberOfRows, numberOfColumns);
         }
 
         /// <summary>
@@ -184,7 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<Complex> CreateVector(int size)
+        public override Vector<Complex> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -1136,204 +1135,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return result;
         }
 
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override Matrix<Complex> Stack(Matrix<Complex> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount);
-            Stack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override void Stack(Matrix<Complex> lower, Matrix<Complex> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != (RowCount + lower.RowCount) || result.ColumnCount != ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < lower.RowCount; i++)
-            {
-                for (var j = 0; j < lower.ColumnCount; j++)
-                {
-                    result.At(i + RowCount, j, lower.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        ///  Concatenates this matrix with the given matrix.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        public override Matrix<Complex> Append(Matrix<Complex> right)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            var result = new SparseMatrix(RowCount, ColumnCount + right.ColumnCount);
-            Append(right, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Concatenates this matrix with the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        public override void Append(Matrix<Complex> right, Matrix<Complex> result)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.ColumnCount != (ColumnCount + right.ColumnCount) || result.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < right.RowCount; i++)
-            {
-                for (var j = 0; j < right.ColumnCount; j++)
-                {
-                    result.At(i, j + RowCount, right.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix, 
-        /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
-        /// The values of off the off diagonal matrices/blocks are set to zero.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <returns>the combined matrix</returns>
-        public override Matrix<Complex> DiagonalStack(Matrix<Complex> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
-            DiagonalStack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<Complex> lower, Matrix<Complex> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != RowCount + lower.RowCount || result.ColumnCount != ColumnCount + lower.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            result.SetSubMatrix(RowCount, lower.RowCount, ColumnCount, lower.ColumnCount, lower);
-        }
-
         /// <summary>
         /// Permute the columns of a matrix according to a permutation.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
index 6e7275a2..50d30bc6 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
@@ -190,16 +190,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <summary>
         /// Creates a <c>SparseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>SparseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new SparseMatrix(numberOfRows, numberOfColumns);
         }
@@ -208,10 +205,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// Creates a <see cref="SparseVector"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="SparseVector"/> with the given dimension.
         /// </returns>
-        public override Vector<Complex> CreateVector(int size)
+        public override Vector<Complex> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -690,56 +688,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
-        /// <summary>
-        /// Diagonally stacks this matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (Rows + lower.rows) x (Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<Complex> lower, Matrix<Complex> result)
-        {
-            var lowerSparseMatrix = lower as SparseMatrix;
-            var resultSparseMatrix = result as SparseMatrix;
-
-            if ((lowerSparseMatrix == null) || (resultSparseMatrix == null))
-            {
-                base.DiagonalStack(lower, result);
-            }
-            else
-            {
-                var resultStorage = resultSparseMatrix.Raw;
-                var lowerStorage = lowerSparseMatrix.Raw;
-
-                if (resultSparseMatrix.RowCount != RowCount + lowerSparseMatrix.RowCount || resultSparseMatrix.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount)
-                {
-                    throw DimensionsDontMatch<ArgumentException>(this, resultSparseMatrix, lowerSparseMatrix);
-                }
-
-                resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount;
-                resultStorage.Values = new Complex[resultStorage.ValueCount];
-                resultStorage.ColumnIndices = new int[resultStorage.ValueCount];
-
-                Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount);
-                Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount);
-
-                Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount);
-                Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount);
-
-                // Copy and adjust lower column indices and rowIndex
-                for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++)
-                {
-                    resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount;
-                }
-
-                for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++)
-                {
-                    resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount;
-                }
-            }
-        }
-
         #region Static constructors for special matrices.
         /// <summary>
         /// Initializes a square <see cref="SparseMatrix"/> with all zero's except for ones on the diagonal.
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 8fe6314f..614f4cd8 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -172,16 +172,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <summary>
         /// Creates a <c>DenseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DenseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new DenseMatrix(numberOfRows, numberOfColumns);
         }
@@ -190,10 +187,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// Creates a <see cref="Vector{T}"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<Complex32> CreateVector(int size)
+        public override Vector<Complex32> CreateVector(int size, bool fullyMutable = false)
         {
             return new DenseVector(size);
         }
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index 400d3f50..4c511431 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
@@ -163,18 +163,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DiagonalMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
-            return new DiagonalMatrix(numberOfRows, numberOfColumns);
+            return fullyMutable
+                ? (Matrix<Complex32>) new SparseMatrix(numberOfRows, numberOfColumns)
+                : new DiagonalMatrix(numberOfRows, numberOfColumns);
         }
 
         /// <summary>
@@ -184,7 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<Complex32> CreateVector(int size)
+        public override Vector<Complex32> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -1136,204 +1135,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return result;
         }
 
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override Matrix<Complex32> Stack(Matrix<Complex32> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount);
-            Stack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override void Stack(Matrix<Complex32> lower, Matrix<Complex32> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != (RowCount + lower.RowCount) || result.ColumnCount != ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, lower, result);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < lower.RowCount; i++)
-            {
-                for (var j = 0; j < lower.ColumnCount; j++)
-                {
-                    result.At(i + RowCount, j, lower.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        ///  Concatenates this matrix with the given matrix.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        public override Matrix<Complex32> Append(Matrix<Complex32> right)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            var result = new SparseMatrix(RowCount, ColumnCount + right.ColumnCount);
-            Append(right, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Concatenates this matrix with the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        public override void Append(Matrix<Complex32> right, Matrix<Complex32> result)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.ColumnCount != (ColumnCount + right.ColumnCount) || result.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < right.RowCount; i++)
-            {
-                for (var j = 0; j < right.ColumnCount; j++)
-                {
-                    result.At(i, j + RowCount, right.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix, 
-        /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
-        /// The values of off the off diagonal matrices/blocks are set to zero.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <returns>the combined matrix</returns>
-        public override Matrix<Complex32> DiagonalStack(Matrix<Complex32> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
-            DiagonalStack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<Complex32> lower, Matrix<Complex32> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != RowCount + lower.RowCount || result.ColumnCount != ColumnCount + lower.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, lower, result);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            result.SetSubMatrix(RowCount, lower.RowCount, ColumnCount, lower.ColumnCount, lower);
-        }
-
         /// <summary>
         /// Permute the columns of a matrix according to a permutation.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
index 8164cb4b..f780d53a 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
@@ -190,16 +190,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <summary>
         /// Creates a <c>SparseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>SparseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new SparseMatrix(numberOfRows, numberOfColumns);
         }
@@ -208,10 +205,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// Creates a <see cref="SparseVector"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="SparseVector"/> with the given dimension.
         /// </returns>
-        public override Vector<Complex32> CreateVector(int size)
+        public override Vector<Complex32> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -690,56 +688,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
-        /// <summary>
-        /// Diagonally stacks this matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (Rows + lower.rows) x (Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<Complex32> lower, Matrix<Complex32> result)
-        {
-            var lowerSparseMatrix = lower as SparseMatrix;
-            var resultSparseMatrix = result as SparseMatrix;
-
-            if ((lowerSparseMatrix == null) || (resultSparseMatrix == null))
-            {
-                base.DiagonalStack(lower, result);
-            }
-            else
-            {
-                var resultStorage = resultSparseMatrix.Raw;
-                var lowerStorage = lowerSparseMatrix.Raw;
-
-                if (resultSparseMatrix.RowCount != RowCount + lowerSparseMatrix.RowCount || resultSparseMatrix.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount)
-                {
-                    throw DimensionsDontMatch<ArgumentException>(this, resultSparseMatrix, lowerSparseMatrix);
-                }
-
-                resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount;
-                resultStorage.Values = new Complex32[resultStorage.ValueCount];
-                resultStorage.ColumnIndices = new int[resultStorage.ValueCount];
-
-                Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount);
-                Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount);
-
-                Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount);
-                Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount);
-
-                // Copy and adjust lower column indices and rowIndex
-                for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++)
-                {
-                    resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount;
-                }
-
-                for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++)
-                {
-                    resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount;
-                }
-            }
-        }
-
         #region Static constructors for special matrices.
         /// <summary>
         /// Initializes a square <see cref="SparseMatrix"/> with all zero's except for ones on the diagonal.
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index 16756af5..50007cd9 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -172,16 +172,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Creates a <c>DenseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DenseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new DenseMatrix(numberOfRows, numberOfColumns);
         }
@@ -190,10 +187,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// Creates a <see cref="Vector{T}"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<double> CreateVector(int size)
+        public override Vector<double> CreateVector(int size, bool fullyMutable = false)
         {
             return new DenseVector(size);
         }
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 4ef3b400..c0da0aa4 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -162,18 +162,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DiagonalMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
-            return new DiagonalMatrix(numberOfRows, numberOfColumns);
+            return fullyMutable
+                ? (Matrix<double>) new SparseMatrix(numberOfRows, numberOfColumns)
+                : new DiagonalMatrix(numberOfRows, numberOfColumns);
         }
 
         /// <summary>
@@ -183,7 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<double> CreateVector(int size)
+        public override Vector<double> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -1130,204 +1129,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return result;
         }
 
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override Matrix<double> Stack(Matrix<double> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount);
-            Stack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override void Stack(Matrix<double> lower, Matrix<double> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != (RowCount + lower.RowCount) || result.ColumnCount != ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, lower, result);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < lower.RowCount; i++)
-            {
-                for (var j = 0; j < lower.ColumnCount; j++)
-                {
-                    result.At(i + RowCount, j, lower.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        ///  Concatenates this matrix with the given matrix.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        public override Matrix<double> Append(Matrix<double> right)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            var result = new SparseMatrix(RowCount, ColumnCount + right.ColumnCount);
-            Append(right, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Concatenates this matrix with the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        public override void Append(Matrix<double> right, Matrix<double> result)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.ColumnCount != (ColumnCount + right.ColumnCount) || result.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < right.RowCount; i++)
-            {
-                for (var j = 0; j < right.ColumnCount; j++)
-                {
-                    result.At(i, j + RowCount, right.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix, 
-        /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
-        /// The values of off the off diagonal matrices/blocks are set to zero.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <returns>the combined matrix</returns>
-        public override Matrix<double> DiagonalStack(Matrix<double> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
-            DiagonalStack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<double> lower, Matrix<double> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != RowCount + lower.RowCount || result.ColumnCount != ColumnCount + lower.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, lower, result);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            result.SetSubMatrix(RowCount, lower.RowCount, ColumnCount, lower.ColumnCount, lower);
-        }
-
         /// <summary>
         /// Permute the columns of a matrix according to a permutation.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index 8b9cf66a..e8cbbc09 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
@@ -189,16 +189,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <summary>
         /// Creates a <c>SparseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>SparseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new SparseMatrix(numberOfRows, numberOfColumns);
         }
@@ -207,10 +204,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// Creates a <see cref="SparseVector"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="SparseVector"/> with the given dimension.
         /// </returns>
-        public override Vector<double> CreateVector(int size)
+        public override Vector<double> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -688,56 +686,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
-        /// <summary>
-        /// Diagonally stacks this matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (Rows + lower.rows) x (Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<double> lower, Matrix<double> result)
-        {
-            var lowerSparseMatrix = lower as SparseMatrix;
-            var resultSparseMatrix = result as SparseMatrix;
-
-            if ((lowerSparseMatrix == null) || (resultSparseMatrix == null))
-            {
-                base.DiagonalStack(lower, result);
-            }
-            else
-            {
-                var resultStorage = resultSparseMatrix.Raw;
-                var lowerStorage = lowerSparseMatrix.Raw;
-
-                if (resultStorage.RowCount != RowCount + lowerStorage.RowCount || resultStorage.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount)
-                {
-                    throw DimensionsDontMatch<ArgumentException>(this, lowerSparseMatrix, resultSparseMatrix);
-                }
-
-                resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount;
-                resultStorage.Values = new double[resultStorage.ValueCount];
-                resultStorage.ColumnIndices = new int[resultStorage.ValueCount];
-
-                Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount);
-                Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount);
-
-                Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount);
-                Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount);
-
-                // Copy and adjust lower column indices and rowIndex
-                for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++)
-                {
-                    resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount;
-                }
-
-                for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++)
-                {
-                    resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount;
-                }
-            }
-        }
-
         #region Static constructors for special matrices.
         /// <summary>
         /// Initializes a square <see cref="SparseMatrix"/> with all zero's except for ones on the diagonal.
diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.cs
index 0c230c1a..d95d2f76 100644
--- a/src/Numerics/LinearAlgebra/Generic/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Generic/Matrix.cs
@@ -261,31 +261,29 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <summary>
         /// Creates a <strong>Matrix</strong> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <strong>Matrix</strong> with the given dimensions.
         /// </returns>
         /// <remarks>
         /// Creates a matrix of the same matrix type as the current matrix.
         /// </remarks>
-        public abstract Matrix<T> CreateMatrix(int numberOfRows, int numberOfColumns);
+        public abstract Matrix<T> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false);
 
         /// <summary>
         /// Creates a Vector with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A Vector with the given dimension.
         /// </returns>
         /// <remarks>
         /// Creates a vector of the same type as the current matrix.
         /// </remarks>
-        public abstract Vector<T> CreateVector(int size);
+        public abstract Vector<T> CreateVector(int size, bool fullyMutable = false);
 
         /// <summary>
         /// Returns a <see cref="System.String"/> that represents this instance.
@@ -1070,7 +1068,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <item>the size of <paramref name="subMatrix"/> is not at least <paramref name="rowCount"/> x <paramref name="columnCount"/>.</item>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
         /// is not positive.</exception>
-        public virtual void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<T> subMatrix)
+        public void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<T> subMatrix)
         {
             if (subMatrix == null)
             {
@@ -1409,7 +1407,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="right">The matrix to concatenate.</param>
         /// <returns>The combined matrix.</returns>
-        public virtual Matrix<T> Append(Matrix<T> right)
+        public Matrix<T> Append(Matrix<T> right)
         {
             if (right == null)
             {
@@ -1421,8 +1419,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
                 throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
             }
 
-            var result = CreateMatrix(RowCount, ColumnCount + right.ColumnCount);
-            Append(right, result);
+            var result = CreateMatrix(RowCount, ColumnCount + right.ColumnCount, fullyMutable: true);
+            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, skipClearing: true);
+            right.Storage.CopySubMatrixTo(result.Storage, 0, 0, right.RowCount, 0, ColumnCount, right.ColumnCount, skipClearing: true);
             return result;
         }
 
@@ -1431,7 +1430,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="right">The matrix to concatenate.</param>
         /// <param name="result">The combined matrix.</param>
-        public virtual void Append(Matrix<T> right, Matrix<T> result)
+        public void Append(Matrix<T> right, Matrix<T> result)
         {
             if (right == null)
             {
@@ -1453,8 +1452,8 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
                 throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
             }
 
-            result.SetSubMatrix(0, RowCount, 0, ColumnCount, this);
-            result.SetSubMatrix(0, right.RowCount, ColumnCount, right.ColumnCount, right);
+            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount);
+            right.Storage.CopySubMatrixTo(result.Storage, 0, 0, right.RowCount, 0, ColumnCount, right.ColumnCount);
         }
 
         /// <summary>
@@ -1464,7 +1463,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <returns>The combined matrix.</returns>
         /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public virtual Matrix<T> Stack(Matrix<T> lower)
+        public Matrix<T> Stack(Matrix<T> lower)
         {
             if (lower == null)
             {
@@ -1476,8 +1475,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
                 throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
             }
 
-            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount);
-            Stack(lower, result);
+            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount, fullyMutable: true);
+            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, skipClearing: true);
+            lower.Storage.CopySubMatrixTo(result.Storage, 0, RowCount, lower.RowCount, 0, 0, lower.ColumnCount, skipClearing: true);
             return result;
         }
 
@@ -1488,7 +1488,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <param name="result">The combined matrix.</param>
         /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public virtual void Stack(Matrix<T> lower, Matrix<T> result)
+        public void Stack(Matrix<T> lower, Matrix<T> result)
         {
             if (lower == null)
             {
@@ -1510,8 +1510,8 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
                 throw DimensionsDontMatch<ArgumentException>(this, result, "result");
             }
 
-            result.SetSubMatrix(0, RowCount, 0, ColumnCount, this);
-            result.SetSubMatrix(RowCount, lower.RowCount, 0, lower.ColumnCount, lower);
+            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount);
+            lower.Storage.CopySubMatrixTo(result.Storage, 0, RowCount, lower.RowCount, 0, 0, lower.ColumnCount);
         }
 
         /// <summary>
@@ -1522,15 +1522,16 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <param name="lower">The lower, right matrix.</param>
         /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
         /// <returns>the combined matrix</returns>
-        public virtual Matrix<T> DiagonalStack(Matrix<T> lower)
+        public Matrix<T> DiagonalStack(Matrix<T> lower)
         {
             if (lower == null)
             {
                 throw new ArgumentNullException("lower");
             }
 
-            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
-            DiagonalStack(lower, result);
+            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount, fullyMutable: true);
+            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount);
+            lower.Storage.CopySubMatrixTo(result.Storage, 0, RowCount, lower.RowCount, 0, ColumnCount, lower.ColumnCount);
             return result;
         }
 
@@ -1542,7 +1543,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
         /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).</exception>
-        public virtual void DiagonalStack(Matrix<T> lower, Matrix<T> result)
+        public void DiagonalStack(Matrix<T> lower, Matrix<T> result)
         {
             if (lower == null)
             {
@@ -1559,8 +1560,8 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
                 throw DimensionsDontMatch<ArgumentException>(this, result, "result");
             }
 
-            result.SetSubMatrix(0, RowCount, 0, ColumnCount, this);
-            result.SetSubMatrix(RowCount, lower.RowCount, ColumnCount, lower.ColumnCount, lower);
+            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount);
+            lower.Storage.CopySubMatrixTo(result.Storage, 0, RowCount, lower.RowCount, 0, ColumnCount, lower.ColumnCount);
         }
 
         /// <summary>Calculates the L1 norm.</summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 2e4c9463..b2b46989 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -172,16 +172,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <summary>
         /// Creates a <c>DenseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DenseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new DenseMatrix(numberOfRows, numberOfColumns);
         }
@@ -190,10 +187,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// Creates a <see cref="Vector{T}"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<float> CreateVector(int size)
+        public override Vector<float> CreateVector(int size, bool fullyMutable = false)
         {
             return new DenseVector(size);
         }
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index aa15c269..0ea96a27 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -162,18 +162,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <summary>
         /// Creates a <c>DiagonalMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>DiagonalMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
-            return new DiagonalMatrix(numberOfRows, numberOfColumns);
+            return fullyMutable
+                ? (Matrix<float>) new SparseMatrix(numberOfRows, numberOfColumns)
+                : new DiagonalMatrix(numberOfRows, numberOfColumns);
         }
 
         /// <summary>
@@ -183,7 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <returns>
         /// A <see cref="Vector{T}"/> with the given dimension.
         /// </returns>
-        public override Vector<float> CreateVector(int size)
+        public override Vector<float> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -1130,204 +1129,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return result;
         }
 
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override Matrix<float> Stack(Matrix<float> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount);
-            Stack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Stacks this matrix on top of the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="lower">The matrix to stack this matrix upon.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If <strong>upper.Columns != lower.Columns</strong>.</exception>
-        public override void Stack(Matrix<float> lower, Matrix<float> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (lower.ColumnCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != (RowCount + lower.RowCount) || result.ColumnCount != ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, lower, result);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < lower.RowCount; i++)
-            {
-                for (var j = 0; j < lower.ColumnCount; j++)
-                {
-                    result.At(i + RowCount, j, lower.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        ///  Concatenates this matrix with the given matrix.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <returns>The combined <see cref="SparseMatrix"/>.</returns>
-        public override Matrix<float> Append(Matrix<float> right)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            var result = new SparseMatrix(RowCount, ColumnCount + right.ColumnCount);
-            Append(right, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Concatenates this matrix with the given matrix and places the result into the result <see cref="SparseMatrix"/>.
-        /// </summary>
-        /// <param name="right">The matrix to concatenate.</param>
-        /// <param name="result">The combined <see cref="SparseMatrix"/>.</param>
-        public override void Append(Matrix<float> right, Matrix<float> result)
-        {
-            if (right == null)
-            {
-                throw new ArgumentNullException("right");
-            }
-
-            if (right.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.ColumnCount != (ColumnCount + right.ColumnCount) || result.RowCount != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            for (var i = 0; i < right.RowCount; i++)
-            {
-                for (var j = 0; j < right.ColumnCount; j++)
-                {
-                    result.At(i, j + RowCount, right.At(i, j));
-                }
-            }
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix, 
-        /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
-        /// The values of off the off diagonal matrices/blocks are set to zero.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <returns>the combined matrix</returns>
-        public override Matrix<float> DiagonalStack(Matrix<float> lower)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            var result = new SparseMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
-            DiagonalStack(lower, result);
-            return result;
-        }
-
-        /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<float> lower, Matrix<float> result)
-        {
-            if (lower == null)
-            {
-                throw new ArgumentNullException("lower");
-            }
-
-            if (result == null)
-            {
-                throw new ArgumentNullException("result");
-            }
-
-            if (result.RowCount != RowCount + lower.RowCount || result.ColumnCount != ColumnCount + lower.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, lower, result);
-            }
-
-            // Clear the result matrix
-            result.Clear();
-
-            // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < _data.Length; i++)
-            {
-                result.At(i, i, _data[i]);
-            }
-
-            // Copy the lower matrix into the result matrix.
-            result.SetSubMatrix(RowCount, lower.RowCount, ColumnCount, lower.ColumnCount, lower);
-        }
-
         /// <summary>
         /// Permute the columns of a matrix according to a permutation.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index 4f1d6b5a..f72424e2 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
@@ -189,16 +189,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <summary>
         /// Creates a <c>SparseMatrix</c> for the given number of rows and columns.
         /// </summary>
-        /// <param name="numberOfRows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="numberOfColumns">
-        /// The number of columns.
-        /// </param>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
         /// <returns>
         /// A <c>SparseMatrix</c> with the given dimensions.
         /// </returns>
-        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new SparseMatrix(numberOfRows, numberOfColumns);
         }
@@ -207,10 +204,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// Creates a <see cref="SparseVector"/> with a the given dimension.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
         /// <returns>
         /// A <see cref="SparseVector"/> with the given dimension.
         /// </returns>
-        public override Vector<float> CreateVector(int size)
+        public override Vector<float> CreateVector(int size, bool fullyMutable = false)
         {
             return new SparseVector(size);
         }
@@ -688,56 +686,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
-        /// <summary>
-        /// Diagonally stacks this matrix on top of the given matrix and places the combined matrix into the result matrix.
-        /// </summary>
-        /// <param name="lower">The lower, right matrix.</param>
-        /// <param name="result">The combined matrix</param>
-        /// <exception cref="ArgumentNullException">If lower is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="ArgumentException">If the result matrix's dimensions are not (Rows + lower.rows) x (Columns + lower.Columns).</exception>
-        public override void DiagonalStack(Matrix<float> lower, Matrix<float> result)
-        {
-            var lowerSparseMatrix = lower as SparseMatrix;
-            var resultSparseMatrix = result as SparseMatrix;
-
-            if ((lowerSparseMatrix == null) || (resultSparseMatrix == null))
-            {
-                base.DiagonalStack(lower, result);
-            }
-            else
-            {
-                var resultStorage = resultSparseMatrix.Raw;
-                var lowerStorage = lowerSparseMatrix.Raw;
-
-                if (resultStorage.RowCount != RowCount + lowerStorage.RowCount || resultStorage.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount)
-                {
-                    throw DimensionsDontMatch<ArgumentException>(this, lowerSparseMatrix, resultSparseMatrix);
-                }
-
-                resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount;
-                resultStorage.Values = new float[resultStorage.ValueCount];
-                resultStorage.ColumnIndices = new int[resultStorage.ValueCount];
-
-                Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount);
-                Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount);
-
-                Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount);
-                Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount);
-
-                // Copy and adjust lower column indices and rowIndex
-                for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++)
-                {
-                    resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount;
-                }
-
-                for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++)
-                {
-                    resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount;
-                }
-            }
-        }
-
         #region Static constructors for special matrices.
         /// <summary>
         /// Initializes a square <see cref="SparseMatrix"/> with all zero's except for ones on the diagonal.
diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
index cdd1f753..f4c9b46a 100644
--- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
@@ -201,7 +201,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Copy(Data, 0, target.Data, 0, Data.Length);
         }
 
-        void CopyTo(SparseCompressedRowMatrixStorage<T> target, bool skipClearing = false)
+        void CopyTo(SparseCompressedRowMatrixStorage<T> target, bool skipClearing)
         {
             if (target == null)
             {
@@ -225,7 +225,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing = false)
+        void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing)
         {
             if (target == null)
             {
@@ -321,7 +321,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         void CopySubMatrixTo(DenseColumnMajorMatrixStorage<T> target,
             int sourceRowIndex, int targetRowIndex, int rowCount,
             int sourceColumnIndex, int targetColumnIndex, int columnCount,
-            bool skipClearing = false)
+            bool skipClearing)
         {
             if (target == null)
             {
@@ -381,7 +381,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         void CopySubMatrixTo(SparseCompressedRowMatrixStorage<T> target,
             int sourceRowIndex, int targetRowIndex, int rowCount,
             int sourceColumnIndex, int targetColumnIndex, int columnCount,
-            bool skipClearing = false)
+            bool skipClearing)
         {
             if (target == null)
             {
diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
index ee0bbea2..a0cc2832 100644
--- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
@@ -1,4 +1,5 @@
 using System;
+using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Storage
@@ -362,11 +363,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 return;
             }
 
-            if (target == null)
-            {
-                throw new ArgumentNullException("target");
-            }
-
             if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
             {
                 var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount);
@@ -385,13 +381,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
-        void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing = false)
+        void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool skipClearing)
         {
-            if (target == null)
-            {
-                throw new ArgumentNullException("target");
-            }
-
             if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
             {
                 var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount);
@@ -427,6 +418,18 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 throw new ArgumentNullException("target");
             }
 
+            var sparseTarget = target as SparseCompressedRowMatrixStorage<T>;
+            if (sparseTarget != null)
+            {
+                CopySubMatrixTo(sparseTarget, 
+                    sourceRowIndex, targetRowIndex, rowCount,
+                    sourceColumnIndex, targetColumnIndex, columnCount,
+                    skipClearing);
+                return;
+            }
+
+            // FALL BACK
+
             if (ReferenceEquals(this, target))
             {
                 throw new NotSupportedException();
@@ -436,8 +439,86 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 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]);
+                    }
+                }
+            }
+        }
+
+        void CopySubMatrixTo(SparseCompressedRowMatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount,
+            bool skipClearing)
+        {
+            if (ReferenceEquals(this, target))
+            {
+                throw new NotSupportedException();
+            }
+
+            ValidateSubMatrixRange(target,
+                sourceRowIndex, targetRowIndex, rowCount,
+                sourceColumnIndex, targetColumnIndex, columnCount);
+
+            var rowOffset = targetRowIndex - sourceRowIndex;
+            var columnOffset = targetColumnIndex - sourceColumnIndex;
+
+            // special case for empty target - much faster
+            if (target.ValueCount == 0)
+            {
+                // note: ValueCount is maximum resulting ValueCount (just using max to avoid internal copying)
+                // resulting arrays will likely be smaller - unless all values fit in the chosen range.
+                var values = new List<T>(ValueCount);
+                var columnIndices = new List<int>(ValueCount);
+                var rowPointers = target.RowPointers;
+
+                for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++)
+                {
+                    rowPointers[i + rowOffset] = values.Count;
+
+                    var startIndex = RowPointers[i];
+                    var endIndex = i < RowPointers.Length - 1 ? RowPointers[i + 1] : ValueCount;
+
+                    // note: we might be able to replace this loop with Array.Copy (perf)
+                    for (int j = startIndex; j < endIndex; j++)
+                    {
+                        // check if the column index is in the range
+                        if ((ColumnIndices[j] >= sourceColumnIndex) && (ColumnIndices[j] < sourceColumnIndex + columnCount))
+                        {
+                            values.Add(Values[j]);
+                            columnIndices.Add(ColumnIndices[j] + columnOffset);
+                        }
+                    }
+                }
+
+                for(int i=targetRowIndex + rowCount; i<rowPointers.Length; i++)
+                {
+                    rowPointers[i] = values.Count;
+                }
+
+                target.ValueCount = values.Count;
+                target.Values = values.ToArray();
+                target.ColumnIndices = columnIndices.ToArray();
+
+                return;
+            }
+
+            // NOTE: potential for more efficient implementation
 
             if (!skipClearing)
             {
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedMatrix.cs b/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedMatrix.cs
index 0803b2f3..de9f71bb 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedMatrix.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedMatrix.cs
@@ -99,7 +99,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
         /// <param name="numberOfRows">The number of rows.</param>
         /// <param name="numberOfColumns">The number of columns.</param>
         /// <returns>A matrix with the given dimensions.</returns>
-        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new UserDefinedMatrix(numberOfRows, numberOfColumns);
         }
@@ -109,7 +109,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
         /// </summary>
         /// <param name="size">The size of the vector.</param>
         /// <returns>A vector with the given dimension.</returns>
-        public override Vector<Complex> CreateVector(int size)
+        public override Vector<Complex> CreateVector(int size, bool fullyMutable = false)
         {
             return new UserDefinedVector(size);
         }
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedMatrix.cs b/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedMatrix.cs
index 6f181705..f3ac48cd 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedMatrix.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedMatrix.cs
@@ -99,7 +99,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
         /// <param name="numberOfRows">The number of rows.</param>
         /// <param name="numberOfColumns">The number of columns.</param>
         /// <returns>A matrix with the given dimensions.</returns>
-        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new UserDefinedMatrix(numberOfRows, numberOfColumns);
         }
@@ -109,7 +109,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
         /// </summary>
         /// <param name="size">The size of the vector.</param>
         /// <returns>A vector with the given dimension.</returns>
-        public override Vector<Complex32> CreateVector(int size)
+        public override Vector<Complex32> CreateVector(int size, bool fullyMutable = false)
         {
             return new UserDefinedVector(size);
         }
diff --git a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrix.cs b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrix.cs
index 902d6aeb..b93b4c43 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrix.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrix.cs
@@ -98,7 +98,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         /// <param name="numberOfRows">The number of rows.</param>
         /// <param name="numberOfColumns">The number of columns.</param>
         /// <returns>A matrix with the given dimensions.</returns>
-        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new UserDefinedMatrix(numberOfRows, numberOfColumns);
         }
@@ -108,7 +108,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         /// </summary>
         /// <param name="size">The size of the vector.</param>
         /// <returns>A vector with the given dimension.</returns>
-        public override Vector<double> CreateVector(int size)
+        public override Vector<double> CreateVector(int size, bool fullyMutable = false)
         {
             return new UserDefinedVector(size);
         }
diff --git a/src/UnitTests/LinearAlgebraTests/Single/UserDefinedMatrix.cs b/src/UnitTests/LinearAlgebraTests/Single/UserDefinedMatrix.cs
index 0679607f..583ae21b 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/UserDefinedMatrix.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/UserDefinedMatrix.cs
@@ -98,7 +98,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
         /// <param name="numberOfRows">The number of rows.</param>
         /// <param name="numberOfColumns">The number of columns.</param>
         /// <returns>A matrix with the given dimensions.</returns>
-        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns)
+        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
         {
             return new UserDefinedMatrix(numberOfRows, numberOfColumns);
         }
@@ -108,7 +108,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
         /// </summary>
         /// <param name="size">The size of the vector.</param>
         /// <returns>A vector with the given dimension.</returns>
-        public override Vector<float> CreateVector(int size)
+        public override Vector<float> CreateVector(int size, bool fullyMutable = false)
         {
             return new UserDefinedVector(size);
         }