diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index 7142e826..43adc606 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -407,6 +407,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             get { return _values; }
         }
 
+        /// <summary>Calculates the induced L1 norm of this matrix.</summary>
+        /// <returns>The maximum absolute column sum of the matrix.</returns>
+        public override double L1Norm()
+        {
+            return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values);
+        }
+
+        /// <summary>Calculates the induced infinity norm of this matrix.</summary>
+        /// <returns>The maximum absolute row sum of the matrix.</returns>
+        public override double InfinityNorm()
+        {
+            return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
+        }
+
+        /// <summary>Calculates the entry-wise Frobenius norm of this matrix.</summary>
+        /// <returns>The square root of the sum of the squared values.</returns>
+        public override double FrobeniusNorm()
+        {
+            return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
+        }
+
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>
@@ -422,29 +443,139 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                     ret._values[(i * _columnCount) + j] = _values[index + i];
                 }
             }
+            return ret;
+        }
 
+        /// <summary>
+        /// Returns the conjugate transpose of this matrix.
+        /// </summary>
+        /// <returns>The conjugate transpose of this matrix.</returns>
+        public override Matrix<Complex> ConjugateTranspose()
+        {
+            var ret = new DenseMatrix(_columnCount, _rowCount);
+            for (var j = 0; j < _columnCount; j++)
+            {
+                var index = j * _rowCount;
+                for (var i = 0; i < _rowCount; i++)
+                {
+                    ret._values[(i * _columnCount) + j] = _values[index + i].Conjugate();
+                }
+            }
             return ret;
         }
 
-        /// <summary>Calculates the induced L1 norm of this matrix.</summary>
-        /// <returns>The maximum absolute column sum of the matrix.</returns>
-        public override double L1Norm()
+        /// <summary>
+        /// Add a scalar to each element of the matrix and stores the result in the result vector.
+        /// </summary>
+        /// <param name="scalar">The scalar to add.</param>
+        /// <param name="result">The matrix to store the result of the addition.</param>
+        protected override void DoAdd(Complex scalar, Matrix<Complex> result)
         {
-            return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values);
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoAdd(scalar, result);
+                return;
+            }
+
+            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] = _values[i] + scalar;
+                }
+            });
         }
 
-        /// <summary>Calculates the induced infinity norm of this matrix.</summary>
-        /// <returns>The maximum absolute row sum of the matrix.</returns>
-        public override double InfinityNorm()
+        /// <summary>
+        /// Adds another matrix to this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to add to this matrix.</param>
+        /// <param name="result">The matrix to store the result of add</param>
+        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
+        protected override void DoAdd(Matrix<Complex> other, Matrix<Complex> result)
         {
-            return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<Complex>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<Complex>;
+            if (denseOther != null && denseResult != null)
+            {
+                Control.LinearAlgebraProvider.AddArrays(_values, denseOther.Data, denseResult.Data);
+                return;
+            }
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex>;
+            if (diagonalOther != null)
+            {
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) + diagonal[i]);
+                }
+                return;
+            }
+
+            base.DoAdd(other, result);
         }
 
-        /// <summary>Calculates the entry-wise Frobenius norm of this matrix.</summary>
-        /// <returns>The square root of the sum of the squared values.</returns>
-        public override double FrobeniusNorm()
+        /// <summary>
+        /// Subtracts a scalar from each element of the matrix and stores the result in the result vector.
+        /// </summary>
+        /// <param name="scalar">The scalar to subtract.</param>
+        /// <param name="result">The matrix to store the result of the subtraction.</param>
+        protected override void DoSubtract(Complex scalar, Matrix<Complex> result)
         {
-            return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoSubtract(scalar, result);
+                return;
+            }
+
+            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] = _values[i] - scalar;
+                }
+            });
+        }
+
+        /// <summary>
+        /// Subtracts another matrix from this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to subtract.</param>
+        /// <param name="result">The matrix to store the result of the subtraction.</param>
+        protected override void DoSubtract(Matrix<Complex> other, Matrix<Complex> result)
+        {
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<Complex>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<Complex>;
+            if (denseOther != null && denseResult != null)
+            {
+                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther.Data, denseResult.Data);
+                return;
+            }
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex>;
+            if (diagonalOther != null)
+            {
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) - diagonal[i]);
+                }
+                return;
+            }
+
+            base.DoSubtract(other, result);
         }
 
         /// <summary>
@@ -505,12 +636,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
             var denseOther = other as DenseMatrix;
             var denseResult = result as DenseMatrix;
-
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoMultiply(other, result);
-            }
-            else
+            if (denseOther != null && denseResult != null)
             {
                 Control.LinearAlgebraProvider.MatrixMultiplyWithUpdate(
                     Providers.LinearAlgebra.Transpose.DontTranspose,
@@ -524,7 +650,31 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                     denseOther._columnCount,
                     0.0,
                     denseResult._values);
+                return;
             }
+
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex>;
+            if (diagonalOther != null)
+            {
+                var diagonal = diagonalOther.Data;
+                var d = Math.Min(ColumnCount, other.ColumnCount);
+                if (d < other.ColumnCount)
+                {
+                    result.ClearSubMatrix(0, RowCount, ColumnCount, other.ColumnCount - ColumnCount);
+                }
+                int index = 0;
+                for (int j = 0; j < d; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        result.At(i, j, _values[index]*diagonal[j]);
+                        index++;
+                    }
+                }
+                return;
+            }
+            
+            base.DoMultiply(other, result);
         }
 
         /// <summary>
@@ -696,113 +846,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
-        /// <summary>
-        /// Add a scalar to each element of the matrix and stores the result in the result vector.
-        /// </summary>
-        /// <param name="scalar">The scalar to add.</param>
-        /// <param name="result">The matrix to store the result of the addition.</param>
-        protected override void DoAdd(Complex scalar, Matrix<Complex> result)
-        {
-            var denseResult = result as DenseMatrix;
-            if (denseResult == null)
-            {
-                base.DoAdd(scalar, result);
-                return;
-            }
-
-            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
-                {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] + scalar;
-                    }
-                });
-        }
-
-        /// <summary>
-        /// Adds another matrix to this matrix.
-        /// </summary>
-        /// <param name="other">The matrix to add to this matrix.</param>
-        /// <param name="result">The matrix to store the result of add</param>
-        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
-        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
-        protected override void DoAdd(Matrix<Complex> other, Matrix<Complex> result)
-        {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoAdd(other, result);
-            }
-            else
-            {
-                Control.LinearAlgebraProvider.AddArrays(_values, denseOther._values, denseResult._values);
-            }
-        }
-
-        /// <summary>
-        /// Subtracts a scalar from each element of the matrix and stores the result in the result vector.
-        /// </summary>
-        /// <param name="scalar">The scalar to subtract.</param>
-        /// <param name="result">The matrix to store the result of the subtraction.</param>
-        protected override void DoSubtract(Complex scalar, Matrix<Complex> result)
-        {
-            var denseResult = result as DenseMatrix;
-            if (denseResult == null)
-            {
-                base.DoSubtract(scalar, result);
-                return;
-            }
-
-            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
-                {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] - scalar;
-                    }
-                });
-        }
-
-        /// <summary>
-        /// Subtracts another matrix from this matrix.
-        /// </summary>
-        /// <param name="other">The matrix to subtract.</param>
-        /// <param name="result">The matrix to store the result of the subtraction.</param>
-        protected override void DoSubtract(Matrix<Complex> other, Matrix<Complex> result)
-        {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoSubtract(other, result);
-            }
-            else
-            {
-                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther._values, denseResult._values);
-            }
-        }
-
-        /// <summary>
-        /// Returns the conjugate transpose of this matrix.
-        /// </summary>
-        /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<Complex> ConjugateTranspose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i].Conjugate();
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index d6a3e968..c89e3e5b 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -588,9 +588,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// i == j (i is the row index, and j is the column index).</remarks>
         public override Vector<Complex> Diagonal()
         {
-            // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
-            // not influence onto diagonal elements
-            return new DenseVector((Complex[])_data.Clone());
+            return new DenseVector(_data).Clone();
         }
 
         /// <summary>
@@ -616,31 +614,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             Multiply(otherDiagonal.Transpose(), result);
         }
 
-        /// <summary>
-        /// Multiplies this matrix with transpose of another matrix and returns the result.
-        /// </summary>
-        /// <param name="other">The matrix to multiply with.</param>
-        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>
-        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
-        /// <returns>The result of multiplication.</returns>
-        public override Matrix<Complex> TransposeAndMultiply(Matrix<Complex> other)
-        {
-            var otherDiagonal = other as DiagonalMatrix;
-            if (otherDiagonal == null)
-            {
-                return base.TransposeAndMultiply(other);
-            }
-
-            if (ColumnCount != otherDiagonal.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, otherDiagonal);
-            }
-
-            var result = other.CreateMatrix(RowCount, other.RowCount);
-            TransposeAndMultiply(other, result);
-            return result;
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>
@@ -703,7 +676,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
             if (RowCount != ColumnCount)
             {
-                    throw new ArgumentException(Resources.ArgumentMatrixSquare);
+                throw new ArgumentException(Resources.ArgumentMatrixSquare);
             }
 
             var inverse = (DiagonalMatrix)Clone();
diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
index 75c1412d..a15cd12c 100644
--- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
@@ -41,6 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
     using Complex = Numerics.Complex;
 #else
     using Complex = System.Numerics.Complex;
+
 #endif
 
     /// <summary>
@@ -48,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
     /// </summary>
     [Serializable]
     public abstract class Matrix : Matrix<Complex>
-    {        
+    {
         /// <summary>
         /// Initializes a new instance of the Matrix class.
         /// </summary>
@@ -96,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public override double FrobeniusNorm()
         {
             var transpose = ConjugateTranspose();
-            var aat = this * transpose;
+            var aat = this*transpose;
             var norm = 0d;
             for (var i = 0; i < RowCount; i++)
             {
@@ -119,7 +120,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                     ret.At(j, i, At(i, j).Conjugate());
                 }
             }
-
             return ret;
         }
 
@@ -202,7 +202,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, At(i, j) * scalar);
+                    result.At(i, j, At(i, j)*scalar);
                 }
             }
         }
@@ -213,18 +213,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <param name="rightSide">The vector to multiply with.</param>
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Vector<Complex> rightSide, Vector<Complex> result)
-         {
+        {
             for (var i = 0; i < RowCount; i++)
             {
                 var s = Complex.Zero;
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
-                    s += At(i, j) * rightSide[j];
+                    s += At(i, j)*rightSide[j];
                 }
-
                 result[i] = s;
             }
-         }
+        }
 
         /// <summary>
         /// Multiplies this matrix with another matrix and places the results into the result matrix.
@@ -240,9 +239,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                     var s = Complex.Zero;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(j, l) * other.At(l, i);
+                        s += At(j, l)*other.At(l, i);
                     }
-
                     result.At(j, i, s);
                 }
             }
@@ -255,7 +253,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <param name="result">The matrix to store the result of the division.</param>
         protected override void DoDivide(Complex divisor, Matrix<Complex> result)
         {
-            DoMultiply(1.0 / divisor, result);
+            DoMultiply(1.0/divisor, result);
         }
 
         /// <summary>
@@ -269,7 +267,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, dividend / At(i, j));
+                    result.At(i, j, dividend/At(i, j));
                 }
             }
         }
@@ -288,9 +286,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                     var s = Complex.Zero;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(i, l) * other.At(j, l);
+                        s += At(i, l)*other.At(j, l);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -310,9 +307,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                     var s = Complex.Zero;
                     for (var l = 0; l < RowCount; l++)
                     {
-                        s += At(l, i) * other.At(l, j);
+                        s += At(l, i)*other.At(l, j);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -328,11 +324,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             for (var i = 0; i < ColumnCount; i++)
             {
                 var s = Complex.Zero;
-                for (var j = 0; j != RowCount; j++)
+                for (var j = 0; j < RowCount; j++)
                 {
-                    s += At(j, i) * rightSide[j];
+                    s += At(j, i)*rightSide[j];
                 }
-
                 result[i] = s;
             }
         }
@@ -345,7 +340,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
             for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     result.At(i, j, -At(i, j));
                 }
@@ -360,7 +355,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
             for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     result.At(i, j, At(i, j).Conjugate());
                 }
@@ -378,7 +373,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             {
                 for (var i = 0; i < RowCount; i++)
                 {
-                    result.At(i, j, At(i, j) * other.At(i, j));
+                    result.At(i, j, At(i, j)*other.At(i, j));
                 }
             }
         }
@@ -394,7 +389,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             {
                 for (var i = 0; i < RowCount; i++)
                 {
-                    result.At(i, j, At(i, j) / divisor.At(i, j));
+                    result.At(i, j, At(i, j)/divisor.At(i, j));
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index df9a2192..37721c1a 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -402,6 +402,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             get { return _values; }
         }
 
+        /// <summary>Calculates the induced L1 norm of this matrix.</summary>
+        /// <returns>The maximum absolute column sum of the matrix.</returns>
+        public override double L1Norm()
+        {
+            return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values);
+        }
+
+        /// <summary>Calculates the induced infinity norm of this matrix.</summary>
+        /// <returns>The maximum absolute row sum of the matrix.</returns>
+        public override double InfinityNorm()
+        {
+            return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
+        }
+
+        /// <summary>Calculates the entry-wise Frobenius norm of this matrix.</summary>
+        /// <returns>The square root of the sum of the squared values.</returns>
+        public override double FrobeniusNorm()
+        {
+            return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
+        }
+
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>
@@ -417,29 +438,139 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                     ret._values[(i * _columnCount) + j] = _values[index + i];
                 }
             }
+            return ret;
+        }
 
+        /// <summary>
+        /// Returns the conjugate transpose of this matrix.
+        /// </summary>
+        /// <returns>The conjugate transpose of this matrix.</returns>
+        public override Matrix<Complex32> ConjugateTranspose()
+        {
+            var ret = new DenseMatrix(_columnCount, _rowCount);
+            for (var j = 0; j < _columnCount; j++)
+            {
+                var index = j * _rowCount;
+                for (var i = 0; i < _rowCount; i++)
+                {
+                    ret._values[(i * _columnCount) + j] = _values[index + i].Conjugate();
+                }
+            }
             return ret;
         }
 
-        /// <summary>Calculates the induced L1 norm of this matrix.</summary>
-        /// <returns>The maximum absolute column sum of the matrix.</returns>
-        public override double L1Norm()
+        /// <summary>
+        /// Add a scalar to each element of the matrix and stores the result in the result vector.
+        /// </summary>
+        /// <param name="scalar">The scalar to add.</param>
+        /// <param name="result">The matrix to store the result of the addition.</param>
+        protected override void DoAdd(Complex32 scalar, Matrix<Complex32> result)
         {
-            return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values);
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoAdd(scalar, result);
+                return;
+            }
+
+            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] = _values[i] + scalar;
+                }
+            });
         }
 
-        /// <summary>Calculates the induced infinity norm of this matrix.</summary>
-        /// <returns>The maximum absolute row sum of the matrix.</returns>
-        public override double InfinityNorm()
+        /// <summary>
+        /// Adds another matrix to this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to add to this matrix.</param>
+        /// <param name="result">The matrix to store the result of add</param>
+        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
+        protected override void DoAdd(Matrix<Complex32> other, Matrix<Complex32> result)
         {
-            return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<Complex32>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<Complex32>;
+            if (denseOther != null && denseResult != null)
+            {
+                Control.LinearAlgebraProvider.AddArrays(_values, denseOther.Data, denseResult.Data);
+                return;
+            }
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex32>;
+            if (diagonalOther != null)
+            {
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) + diagonal[i]);
+                }
+                return;
+            }
+
+            base.DoAdd(other, result);
         }
 
-        /// <summary>Calculates the entry-wise Frobenius norm of this matrix.</summary>
-        /// <returns>The square root of the sum of the squared values.</returns>
-        public override double FrobeniusNorm()
+        /// <summary>
+        /// Subtracts a scalar from each element of the matrix and stores the result in the result vector.
+        /// </summary>
+        /// <param name="scalar">The scalar to subtract.</param>
+        /// <param name="result">The matrix to store the result of the subtraction.</param>
+        protected override void DoSubtract(Complex32 scalar, Matrix<Complex32> result)
         {
-            return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoSubtract(scalar, result);
+                return;
+            }
+
+            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] = _values[i] - scalar;
+                }
+            });
+        }
+
+        /// <summary>
+        /// Subtracts another matrix from this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to subtract.</param>
+        /// <param name="result">The matrix to store the result of the subtraction.</param>
+        protected override void DoSubtract(Matrix<Complex32> other, Matrix<Complex32> result)
+        {
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<Complex32>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<Complex32>;
+            if (denseOther != null && denseResult != null)
+            {
+                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther.Data, denseResult.Data);
+                return;
+            }
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex32>;
+            if (diagonalOther != null)
+            {
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) - diagonal[i]);
+                }
+                return;
+            }
+
+            base.DoSubtract(other, result);
         }
 
         /// <summary>
@@ -500,12 +631,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         {
             var denseOther = other as DenseMatrix;
             var denseResult = result as DenseMatrix;
-
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoMultiply(other, result);
-            }
-            else
+            if (denseOther != null && denseResult != null)
             {
                 Control.LinearAlgebraProvider.MatrixMultiplyWithUpdate(
                     Providers.LinearAlgebra.Transpose.DontTranspose,
@@ -519,7 +645,31 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                     denseOther._columnCount,
                     0.0f,
                     denseResult._values);
+                return;
             }
+
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex32>;
+            if (diagonalOther != null)
+            {
+                var diagonal = diagonalOther.Data;
+                var d = Math.Min(ColumnCount, other.ColumnCount);
+                if (d < other.ColumnCount)
+                {
+                    result.ClearSubMatrix(0, RowCount, ColumnCount, other.ColumnCount - ColumnCount);
+                }
+                int index = 0;
+                for (int j = 0; j < d; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        result.At(i, j, _values[index]*diagonal[j]);
+                        index++;
+                    }
+                }
+                return;
+            }
+
+            base.DoMultiply(other, result);
         }
 
         /// <summary>
@@ -691,113 +841,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
-        /// <summary>
-        /// Add a scalar to each element of the matrix and stores the result in the result vector.
-        /// </summary>
-        /// <param name="scalar">The scalar to add.</param>
-        /// <param name="result">The matrix to store the result of the addition.</param>
-        protected override void DoAdd(Complex32 scalar, Matrix<Complex32> result)
-        {
-            var denseResult = result as DenseMatrix;
-            if (denseResult == null)
-            {
-                base.DoAdd(scalar, result);
-                return;
-            }
-
-            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
-                {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] + scalar;
-                    }
-                });
-        }
-
-        /// <summary>
-        /// Adds another matrix to this matrix.
-        /// </summary>
-        /// <param name="other">The matrix to add to this matrix.</param>
-        /// <param name="result">The matrix to store the result of add</param>
-        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
-        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
-        protected override void DoAdd(Matrix<Complex32> other, Matrix<Complex32> result)
-        {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoAdd(other, result);
-            }
-            else
-            {
-                Control.LinearAlgebraProvider.AddArrays(_values, denseOther._values, denseResult._values);
-            }
-        }
-
-        /// <summary>
-        /// Subtracts a scalar from each element of the matrix and stores the result in the result vector.
-        /// </summary>
-        /// <param name="scalar">The scalar to subtract.</param>
-        /// <param name="result">The matrix to store the result of the subtraction.</param>
-        protected override void DoSubtract(Complex32 scalar, Matrix<Complex32> result)
-        {
-            var denseResult = result as DenseMatrix;
-            if (denseResult == null)
-            {
-                base.DoSubtract(scalar, result);
-                return;
-            }
-
-            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
-                {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] - scalar;
-                    }
-                });
-        }
-
-        /// <summary>
-        /// Subtracts another matrix from this matrix.
-        /// </summary>
-        /// <param name="other">The matrix to subtract.</param>
-        /// <param name="result">The matrix to store the result of the subtraction.</param>
-        protected override void DoSubtract(Matrix<Complex32> other, Matrix<Complex32> result)
-        {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoSubtract(other, result);
-            }
-            else
-            {
-                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther._values, denseResult._values);
-            }
-        }
-
-        /// <summary>
-        /// Returns the conjugate transpose of this matrix.
-        /// </summary>
-        /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<Complex32> ConjugateTranspose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i].Conjugate();
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index f1f04d3b..2213a10e 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
@@ -583,9 +583,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// i == j (i is the row index, and j is the column index).</remarks>
         public override Vector<Complex32> Diagonal()
         {
-            // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
-            // not influence onto diagonal elements
-            return new DenseVector((Complex32[])_data.Clone());
+            return new DenseVector(_data).Clone();
         }
 
         /// <summary>
@@ -611,31 +609,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             Multiply(otherDiagonal.Transpose(), result);
         }
 
-        /// <summary>
-        /// Multiplies this matrix with transpose of another matrix and returns the result.
-        /// </summary>
-        /// <param name="other">The matrix to multiply with.</param>
-        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>
-        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
-        /// <returns>The result of multiplication.</returns>
-        public override Matrix<Complex32> TransposeAndMultiply(Matrix<Complex32> other)
-        {
-            var otherDiagonal = other as DiagonalMatrix;
-            if (otherDiagonal == null)
-            {
-                return base.TransposeAndMultiply(other);
-            }
-
-            if (ColumnCount != otherDiagonal.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, otherDiagonal);
-            }
-
-            var result = other.CreateMatrix(RowCount, other.RowCount);
-            TransposeAndMultiply(other, result);
-            return result;
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
index ad6a3346..dae6f558 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
@@ -28,11 +28,11 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System;
 using MathNet.Numerics.LinearAlgebra.Complex32.Factorization;
 using MathNet.Numerics.LinearAlgebra.Factorization;
 using MathNet.Numerics.LinearAlgebra.Storage;
 using MathNet.Numerics.Properties;
-using System;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32
 {
@@ -91,7 +91,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public override double FrobeniusNorm()
         {
             var transpose = ConjugateTranspose();
-            var aat = this * transpose;
+            var aat = this*transpose;
             var norm = 0d;
             for (var i = 0; i < RowCount; i++)
             {
@@ -196,29 +196,28 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, At(i, j) * scalar);
+                    result.At(i, j, At(i, j)*scalar);
                 }
             }
         }
 
-         /// <summary>
+        /// <summary>
         /// Multiplies this matrix with a vector and places the results into the result vector.
         /// </summary>
         /// <param name="rightSide">The vector to multiply with.</param>
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Vector<Complex32> rightSide, Vector<Complex32> result)
-         {
+        {
             for (var i = 0; i < RowCount; i++)
             {
                 var s = Complex32.Zero;
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
-                    s += At(i, j) * rightSide[j];
+                    s += At(i, j)*rightSide[j];
                 }
-
                 result[i] = s;
             }
-         }
+        }
 
         /// <summary>
         /// Divides each element of the matrix by a scalar and places results into the result matrix.
@@ -227,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <param name="result">The matrix to store the result of the division.</param>
         protected override void DoDivide(Complex32 divisor, Matrix<Complex32> result)
         {
-            DoMultiply(1.0f / divisor, result);
+            DoMultiply(1.0f/divisor, result);
         }
 
         /// <summary>
@@ -241,7 +240,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, dividend / At(i, j));
+                    result.At(i, j, dividend/At(i, j));
                 }
             }
         }
@@ -260,9 +259,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                     var s = Complex32.Zero;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(j, l) * other.At(l, i);
+                        s += At(j, l)*other.At(l, i);
                     }
-
                     result.At(j, i, s);
                 }
             }
@@ -282,9 +280,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                     var s = Complex32.Zero;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(i, l) * other.At(j, l);
+                        s += At(i, l)*other.At(j, l);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -304,9 +301,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                     var s = Complex32.Zero;
                     for (var l = 0; l < RowCount; l++)
                     {
-                        s += At(l, i) * other.At(l, j);
+                        s += At(l, i)*other.At(l, j);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -322,11 +318,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             for (var i = 0; i < ColumnCount; i++)
             {
                 var s = Complex32.Zero;
-                for (var j = 0; j != RowCount; j++)
+                for (var j = 0; j < RowCount; j++)
                 {
-                    s += At(j, i) * rightSide[j];
+                    s += At(j, i)*rightSide[j];
                 }
-
                 result[i] = s;
             }
         }
@@ -339,7 +334,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         {
             for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     result.At(i, j, -At(i, j));
                 }
@@ -354,7 +349,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         {
             for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     result.At(i, j, At(i, j).Conjugate());
                 }
@@ -372,7 +367,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             {
                 for (var i = 0; i < RowCount; i++)
                 {
-                    result.At(i, j, At(i, j) * other.At(i, j));
+                    result.At(i, j, At(i, j)*other.At(i, j));
                 }
             }
         }
@@ -388,7 +383,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             {
                 for (var i = 0; i < RowCount; i++)
                 {
-                    result.At(i, j, At(i, j) / divisor.At(i, j));
+                    result.At(i, j, At(i, j)/divisor.At(i, j));
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index 30198a63..f9370001 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -399,25 +399,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             get { return _values; }
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<double> Transpose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i];
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()
@@ -439,6 +420,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
         }
 
+        /// <summary>
+        /// Returns the transpose of this matrix.
+        /// </summary>
+        /// <returns>The transpose of this matrix.</returns>
+        public override Matrix<double> Transpose()
+        {
+            var ret = new DenseMatrix(_columnCount, _rowCount);
+            for (var j = 0; j < _columnCount; j++)
+            {
+                var index = j * _rowCount;
+                for (var i = 0; i < _rowCount; i++)
+                {
+                    ret._values[(i * _columnCount) + j] = _values[index + i];
+                }
+            }
+
+            return ret;
+        }
+
         /// <summary>
         /// Add a scalar to each element of the matrix and stores the result in the result vector.
         /// </summary>
@@ -454,13 +454,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] + scalar;
-                    }
-                });
+                    v[i] = _values[i] + scalar;
+                }
+            });
         }
 
         /// <summary>
@@ -472,16 +472,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoAdd(Matrix<double> other, Matrix<double> result)
         {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<double>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<double>;
+            if (denseOther != null && denseResult != null)
             {
-                base.DoAdd(other, result);
+                Control.LinearAlgebraProvider.AddArrays(_values, denseOther.Data, denseResult.Data);
+                return;
             }
-            else
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<double>;
+            if (diagonalOther != null)
             {
-                Control.LinearAlgebraProvider.AddArrays(_values, denseOther._values, denseResult._values);
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) + diagonal[i]);
+                }
+                return;
             }
+
+            base.DoAdd(other, result);
         }
 
         /// <summary>
@@ -499,13 +512,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] - scalar;
-                    }
-                });
+                    v[i] = _values[i] - scalar;
+                }
+            });
         }
 
         /// <summary>
@@ -515,16 +528,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="result">The matrix to store the result of the subtraction.</param>
         protected override void DoSubtract(Matrix<double> other, Matrix<double> result)
         {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<double>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<double>;
+            if (denseOther != null && denseResult != null)
             {
-                base.DoSubtract(other, result);
+                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther.Data, denseResult.Data);
+                return;
             }
-            else
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<double>;
+            if (diagonalOther != null)
             {
-                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther._values, denseResult._values);
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) - diagonal[i]);
+                }
+                return;
             }
+
+            base.DoSubtract(other, result);
         }
 
         /// <summary>
@@ -585,12 +611,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             var denseOther = other as DenseMatrix;
             var denseResult = result as DenseMatrix;
-
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoMultiply(other, result);
-            }
-            else
+            if (denseOther != null && denseResult != null)
             {
                 Control.LinearAlgebraProvider.MatrixMultiplyWithUpdate(
                     Providers.LinearAlgebra.Transpose.DontTranspose,
@@ -604,7 +625,31 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                     denseOther._columnCount,
                     0.0,
                     denseResult._values);
+                return;
             }
+
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<double>;
+            if (diagonalOther != null)
+            {
+                var diagonal = diagonalOther.Data;
+                var d = Math.Min(ColumnCount, other.ColumnCount);
+                if (d < other.ColumnCount)
+                {
+                    result.ClearSubMatrix(0, RowCount, ColumnCount, other.ColumnCount - ColumnCount);
+                }
+                int index = 0;
+                for (int j = 0; j < d; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        result.At(i, j, _values[index]*diagonal[j]);
+                        index++;
+                    }
+                }
+                return;
+            }
+
+            base.DoMultiply(other, result);
         }
 
         /// <summary>
@@ -614,9 +659,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoTransposeAndMultiply(Matrix<double> other, Matrix<double> result)
         {
-             var denseOther = other as DenseMatrix;
+            var denseOther = other as DenseMatrix;
             var denseResult = result as DenseMatrix;
-
             if (denseOther == null || denseResult == null)
             {
                 base.DoTransposeAndMultiply(other, result);
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index e489a675..872d85b3 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -577,9 +577,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// i == j (i is the row index, and j is the column index).</remarks>
         public override Vector<double> Diagonal()
         {
-            // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
-            // not influence onto diagonal elements
-            return new DenseVector((double[])_data.Clone());
+            return new DenseVector(_data).Clone();
         }
 
         /// <summary>
@@ -605,31 +603,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             Multiply(otherDiagonal.Transpose(), result);
         }
 
-        /// <summary>
-        /// Multiplies this matrix with transpose of another matrix and returns the result.
-        /// </summary>
-        /// <param name="other">The matrix to multiply with.</param>
-        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>
-        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
-        /// <returns>The result of multiplication.</returns>
-        public override Matrix<double> TransposeAndMultiply(Matrix<double> other)
-        {
-            var otherDiagonal = other as DiagonalMatrix;
-            if (otherDiagonal == null)
-            {
-                return base.TransposeAndMultiply(other);
-            }
-
-            if (ColumnCount != otherDiagonal.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, otherDiagonal);
-            }
-
-            var result = other.CreateMatrix(RowCount, other.RowCount);
-            TransposeAndMultiply(other, result);
-            return result;
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>
@@ -692,7 +665,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             if (RowCount != ColumnCount)
             {
-                    throw new ArgumentException(Resources.ArgumentMatrixSquare);
+                throw new ArgumentException(Resources.ArgumentMatrixSquare);
             }
 
             var inverse = (DiagonalMatrix)Clone();
diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs
index ec4eaa25..b70d67f7 100644
--- a/src/Numerics/LinearAlgebra/Double/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs
@@ -41,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
     /// </summary>
     [Serializable]
     public abstract class Matrix : Matrix<double>
-    {        
+    {
         /// <summary>
         /// Initializes a new instance of the Matrix class.
         /// </summary>
@@ -89,7 +89,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public override double FrobeniusNorm()
         {
             var transpose = Transpose();
-            var aat = this * transpose;
+            var aat = this*transpose;
             var norm = 0d;
             for (var i = 0; i < RowCount; i++)
             {
@@ -156,7 +156,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 }
             }
         }
-    
+
         /// <summary>
         /// Subtracts another matrix from this matrix.
         /// </summary>
@@ -186,7 +186,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, At(i, j) * scalar);
+                    result.At(i, j, At(i, j)*scalar);
                 }
             }
         }
@@ -197,18 +197,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="rightSide">The vector to multiply with.</param>
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Vector<double> rightSide, Vector<double> result)
-         {
-             for (var i = 0; i < RowCount; i++)
-             {
-                 var s = 0.0;
-                 for (var j = 0; j != ColumnCount; j++)
-                 {
-                     s += At(i, j) * rightSide[j];
-                 }
-
-                 result[i] = s;
-             }
-         }
+        {
+            for (var i = 0; i < RowCount; i++)
+            {
+                var s = 0.0;
+                for (var j = 0; j < ColumnCount; j++)
+                {
+                    s += At(i, j)*rightSide[j];
+                }
+                result[i] = s;
+            }
+        }
 
         /// <summary>
         /// Divides each element of the matrix by a scalar and places results into the result matrix.
@@ -217,7 +216,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="result">The matrix to store the result of the division.</param>
         protected override void DoDivide(double divisor, Matrix<double> result)
         {
-            DoMultiply(1.0 / divisor, result);
+            DoMultiply(1.0/divisor, result);
         }
 
         /// <summary>
@@ -231,7 +230,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, dividend / At(i, j));
+                    result.At(i, j, dividend/At(i, j));
                 }
             }
         }
@@ -245,14 +244,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             for (var j = 0; j < RowCount; j++)
             {
-                for (var i = 0; i != other.ColumnCount; i++)
+                for (var i = 0; i < other.ColumnCount; i++)
                 {
                     var s = 0.0;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(j, l) * other.At(l, i);
+                        s += At(j, l)*other.At(l, i);
                     }
-
                     result.At(j, i, s);
                 }
             }
@@ -272,9 +270,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                     var s = 0.0;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(i, l) * other.At(j, l);
+                        s += At(i, l)*other.At(j, l);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -294,9 +291,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                     var s = 0.0;
                     for (var l = 0; l < RowCount; l++)
                     {
-                        s += At(l, i) * other.At(l, j);
+                        s += At(l, i)*other.At(l, j);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -312,11 +308,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             for (var i = 0; i < ColumnCount; i++)
             {
                 var s = 0.0;
-                for (var j = 0; j != RowCount; j++)
+                for (var j = 0; j < RowCount; j++)
                 {
-                    s += At(j, i) * rightSide[j];
+                    s += At(j, i)*rightSide[j];
                 }
-
                 result[i] = s;
             }
         }
@@ -329,7 +324,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     result.At(i, j, -At(i, j));
                 }
diff --git a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs
index 25163d6a..bb5fcebf 100644
--- a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs
+++ b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs
@@ -702,7 +702,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <param name="other">The matrix to multiply with.</param>
         /// <exception cref="ArgumentException">If <strong>this.Columns != other.ColumnCount</strong>.</exception>
         /// <returns>The result of the multiplication.</returns>
-        public virtual Matrix<T> TransposeAndMultiply(Matrix<T> other)
+        public Matrix<T> TransposeAndMultiply(Matrix<T> other)
         {
             if (ColumnCount != other.ColumnCount)
             {
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 8c3417ee..b112612d 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -399,25 +399,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             get { return _values; }
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<float> Transpose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i];
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()
@@ -439,6 +420,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
         }
 
+        /// <summary>
+        /// Returns the transpose of this matrix.
+        /// </summary>
+        /// <returns>The transpose of this matrix.</returns>
+        public override Matrix<float> Transpose()
+        {
+            var ret = new DenseMatrix(_columnCount, _rowCount);
+            for (var j = 0; j < _columnCount; j++)
+            {
+                var index = j * _rowCount;
+                for (var i = 0; i < _rowCount; i++)
+                {
+                    ret._values[(i * _columnCount) + j] = _values[index + i];
+                }
+            }
+
+            return ret;
+        }
+
         /// <summary>
         /// Add a scalar to each element of the matrix and stores the result in the result vector.
         /// </summary>
@@ -454,13 +454,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] + scalar;
-                    }
-                });
+                    v[i] = _values[i] + scalar;
+                }
+            });
         }
 
         /// <summary>
@@ -472,16 +472,29 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoAdd(Matrix<float> other, Matrix<float> result)
         {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<float>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<float>;
+            if (denseOther != null && denseResult != null)
             {
-                base.DoAdd(other, result);
+                Control.LinearAlgebraProvider.AddArrays(_values, denseOther.Data, denseResult.Data);
+                return;
             }
-            else
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<float>;
+            if (diagonalOther != null)
             {
-                Control.LinearAlgebraProvider.AddArrays(_values, denseOther._values, denseResult._values);
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) + diagonal[i]);
+                }
+                return;
             }
+
+            base.DoAdd(other, result);
         }
 
         /// <summary>
@@ -499,13 +512,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = _values[i] - scalar;
-                    }
-                });
+                    v[i] = _values[i] - scalar;
+                }
+            });
         }
 
         /// <summary>
@@ -515,16 +528,29 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <param name="result">The matrix to store the result of the subtraction.</param>
         protected override void DoSubtract(Matrix<float> other, Matrix<float> result)
         {
-            var denseOther = other as DenseMatrix;
-            var denseResult = result as DenseMatrix;
-            if (denseOther == null || denseResult == null)
+            // dense + dense = dense
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<float>;
+            var denseResult = result.Storage as DenseColumnMajorMatrixStorage<float>;
+            if (denseOther != null && denseResult != null)
             {
-                base.DoSubtract(other, result);
+                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther.Data, denseResult.Data);
+                return;
             }
-            else
+
+            // dense + diagonal = matrix
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<float>;
+            if (diagonalOther != null)
             {
-                Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther._values, denseResult._values);
+                CopyTo(result);
+                var diagonal = diagonalOther.Data;
+                for (int i = 0; i < diagonal.Length; i++)
+                {
+                    result.At(i, i, result.At(i, i) - diagonal[i]);
+                }
+                return;
             }
+
+            base.DoSubtract(other, result);
         }
 
         /// <summary>
@@ -585,12 +611,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
             var denseOther = other as DenseMatrix;
             var denseResult = result as DenseMatrix;
-
-            if (denseOther == null || denseResult == null)
-            {
-                base.DoMultiply(other, result);
-            }
-            else
+            if (denseOther != null && denseResult != null)
             {
                 Control.LinearAlgebraProvider.MatrixMultiplyWithUpdate(
                     Providers.LinearAlgebra.Transpose.DontTranspose,
@@ -604,7 +625,31 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                     denseOther._columnCount,
                     0.0f,
                     denseResult._values);
+                return;
             }
+
+            var diagonalOther = other.Storage as DiagonalMatrixStorage<float>;
+            if (diagonalOther != null)
+            {
+                var diagonal = diagonalOther.Data;
+                var d = Math.Min(ColumnCount, other.ColumnCount);
+                if (d < other.ColumnCount)
+                {
+                    result.ClearSubMatrix(0, RowCount, ColumnCount, other.ColumnCount - ColumnCount);
+                }
+                int index = 0;
+                for (int j = 0; j < d; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        result.At(i, j, _values[index]*diagonal[j]);
+                        index++;
+                    }
+                }
+                return;
+            }
+
+            base.DoMultiply(other, result);
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index e388e2b2..42d5b944 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -577,9 +577,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// i == j (i is the row index, and j is the column index).</remarks>
         public override Vector<float> Diagonal()
         {
-            // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
-            // not influence onto diagonal elements
-            return new DenseVector((float[])_data.Clone());
+            return new DenseVector(_data).Clone();
         }
 
         /// <summary>
@@ -605,31 +603,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             Multiply(otherDiagonal.Transpose(), result);
         }
 
-        /// <summary>
-        /// Multiplies this matrix with transpose of another matrix and returns the result.
-        /// </summary>
-        /// <param name="other">The matrix to multiply with.</param>
-        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>
-        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
-        /// <returns>The result of multiplication.</returns>
-        public override Matrix<float> TransposeAndMultiply(Matrix<float> other)
-        {
-            var otherDiagonal = other as DiagonalMatrix;
-            if (otherDiagonal == null)
-            {
-                return base.TransposeAndMultiply(other);
-            }
-
-            if (ColumnCount != otherDiagonal.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, otherDiagonal);
-            }
-
-            var result = other.CreateMatrix(RowCount, other.RowCount);
-            TransposeAndMultiply(other, result);
-            return result;
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs
index 761734f0..b8737c0e 100644
--- a/src/Numerics/LinearAlgebra/Single/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs
@@ -41,7 +41,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
     /// </summary>
     [Serializable]
     public abstract class Matrix : Matrix<float>
-    {        
+    {
         /// <summary>
         /// Initializes a new instance of the Matrix class.
         /// </summary>
@@ -89,7 +89,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public override double FrobeniusNorm()
         {
             var transpose = Transpose();
-            var aat = this * transpose;
+            var aat = this*transpose;
             var norm = 0d;
             for (var i = 0; i < RowCount; i++)
             {
@@ -186,29 +186,28 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, At(i, j) * scalar);
+                    result.At(i, j, At(i, j)*scalar);
                 }
             }
         }
 
-         /// <summary>
+        /// <summary>
         /// Multiplies this matrix with a vector and places the results into the result vector.
         /// </summary>
         /// <param name="rightSide">The vector to multiply with.</param>
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Vector<float> rightSide, Vector<float> result)
-         {
+        {
             for (var i = 0; i < RowCount; i++)
             {
                 var s = 0.0f;
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
-                    s += At(i, j) * rightSide[j];
+                    s += At(i, j)*rightSide[j];
                 }
-
                 result[i] = s;
             }
-         }
+        }
 
         /// <summary>
         /// Multiplies this matrix with another matrix and places the results into the result matrix.
@@ -219,14 +218,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
             for (var j = 0; j < RowCount; j++)
             {
-                for (var i = 0; i != other.ColumnCount; i++)
+                for (var i = 0; i < other.ColumnCount; i++)
                 {
                     var s = 0.0f;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(j, l) * other.At(l, i);
+                        s += At(j, l)*other.At(l, i);
                     }
-
                     result.At(j, i, s);
                 }
             }
@@ -239,7 +237,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <param name="result">The matrix to store the result of the division.</param>
         protected override void DoDivide(float divisor, Matrix<float> result)
         {
-            DoMultiply(1.0f / divisor, result);
+            DoMultiply(1.0f/divisor, result);
         }
 
         /// <summary>
@@ -253,7 +251,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    result.At(i, j, dividend / At(i, j));
+                    result.At(i, j, dividend/At(i, j));
                 }
             }
         }
@@ -272,9 +270,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                     var s = 0.0f;
                     for (var l = 0; l < ColumnCount; l++)
                     {
-                        s += At(i, l) * other.At(j, l);
+                        s += At(i, l)*other.At(j, l);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -294,9 +291,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                     var s = 0.0f;
                     for (var l = 0; l < RowCount; l++)
                     {
-                        s += At(l, i) * other.At(l, j);
+                        s += At(l, i)*other.At(l, j);
                     }
-
                     result.At(i, j, s);
                 }
             }
@@ -312,11 +308,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             for (var i = 0; i < ColumnCount; i++)
             {
                 var s = 0.0f;
-                for (var j = 0; j != RowCount; j++)
+                for (var j = 0; j < RowCount; j++)
                 {
-                    s += At(j, i) * rightSide[j];
+                    s += At(j, i)*rightSide[j];
                 }
-
                 result[i] = s;
             }
         }
@@ -332,7 +327,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, At(row, column) % divisor);
+                    result.At(row, column, At(row, column)%divisor);
                 }
             }
         }
@@ -348,7 +343,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, dividend % At(row, column));
+                    result.At(row, column, dividend%At(row, column));
                 }
             }
         }
@@ -361,7 +356,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
             for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j != ColumnCount; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     result.At(i, j, -At(i, j));
                 }
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/DiagonalMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/DiagonalMatrixTests.cs
index fc265c89..c54dc08c 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/DiagonalMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/DiagonalMatrixTests.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using MathNet.Numerics.Distributions;
 using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.LinearAlgebra.Complex;
+using MathNet.Numerics.Random;
 using NUnit.Framework;
 
 namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
@@ -383,5 +385,23 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
             var matrix = TestMatrices["Square3x3"];
             Assert.IsTrue(matrix.IsSymmetric);
         }
+
+        [Test]
+        public void DenseDiagonalMatrixMultiplication()
+        {
+            var dist = new ContinuousUniform(-1.0, 1.0, new MersenneTwister());
+
+            Assert.IsInstanceOf<DiagonalMatrix>(Matrix<Complex>.Build.DiagonalIdentity(3, 3));
+
+            var tall = Matrix<Complex>.Build.Random(8, 3, dist);
+            Assert.IsTrue((tall*Matrix<Complex>.Build.DiagonalIdentity(3).Multiply(2d)).Equals(tall.Multiply(2d)));
+            Assert.IsTrue((tall*Matrix<Complex>.Build.Diagonal(3, 5, 2d)).Equals(tall.Multiply(2d).Append(Matrix<Complex>.Build.Dense(8, 2))));
+            Assert.IsTrue((tall*Matrix<Complex>.Build.Diagonal(3, 2, 2d)).Equals(tall.Multiply(2d).SubMatrix(0, 8, 0, 2)));
+
+            var wide = Matrix<Complex>.Build.Random(3, 8, dist);
+            Assert.IsTrue((wide*Matrix<Complex>.Build.DiagonalIdentity(8).Multiply(2d)).Equals(wide.Multiply(2d)));
+            Assert.IsTrue((wide*Matrix<Complex>.Build.Diagonal(8, 10, 2d)).Equals(wide.Multiply(2d).Append(Matrix<Complex>.Build.Dense(3, 2))));
+            Assert.IsTrue((wide*Matrix<Complex>.Build.Diagonal(8, 2, 2d)).Equals(wide.Multiply(2d).SubMatrix(0, 3, 0, 2)));
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/DiagonalMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/DiagonalMatrixTests.cs
index 1e3058e3..2623229d 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/DiagonalMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/DiagonalMatrixTests.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using MathNet.Numerics.Distributions;
 using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.LinearAlgebra.Complex32;
+using MathNet.Numerics.Random;
 using NUnit.Framework;
 
 namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
@@ -379,5 +381,22 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
             var matrix = TestMatrices["Square3x3"];
             Assert.IsTrue(matrix.IsSymmetric);
         }
+        [Test]
+        public void DenseDiagonalMatrixMultiplication()
+        {
+            var dist = new ContinuousUniform(-1.0, 1.0, new MersenneTwister());
+
+            Assert.IsInstanceOf<DiagonalMatrix>(Matrix<Complex32>.Build.DiagonalIdentity(3, 3));
+
+            var tall = Matrix<Complex32>.Build.Random(8, 3, dist);
+            Assert.IsTrue((tall*Matrix<Complex32>.Build.DiagonalIdentity(3).Multiply(2f)).Equals(tall.Multiply(2f)));
+            Assert.IsTrue((tall*Matrix<Complex32>.Build.Diagonal(3, 5, 2f)).Equals(tall.Multiply(2f).Append(Matrix<Complex32>.Build.Dense(8, 2))));
+            Assert.IsTrue((tall*Matrix<Complex32>.Build.Diagonal(3, 2, 2f)).Equals(tall.Multiply(2f).SubMatrix(0, 8, 0, 2)));
+
+            var wide = Matrix<Complex32>.Build.Random(3, 8, dist);
+            Assert.IsTrue((wide*Matrix<Complex32>.Build.DiagonalIdentity(8).Multiply(2f)).Equals(wide.Multiply(2f)));
+            Assert.IsTrue((wide*Matrix<Complex32>.Build.Diagonal(8, 10, 2f)).Equals(wide.Multiply(2f).Append(Matrix<Complex32>.Build.Dense(3, 2))));
+            Assert.IsTrue((wide*Matrix<Complex32>.Build.Diagonal(8, 2, 2f)).Equals(wide.Multiply(2f).SubMatrix(0, 3, 0, 2)));
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Double/DiagonalMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/DiagonalMatrixTests.cs
index 171cfe35..cc3e5240 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/DiagonalMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/DiagonalMatrixTests.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using MathNet.Numerics.Distributions;
 using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Random;
 using NUnit.Framework;
 
 namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
@@ -410,5 +412,23 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             var test = diagonal*dense;
             var test2 = dense*diagonal;
         }
+
+        [Test]
+        public void DenseDiagonalMatrixMultiplication()
+        {
+            var dist = new ContinuousUniform(-1.0, 1.0, new MersenneTwister());
+
+            Assert.IsInstanceOf<DiagonalMatrix>(Matrix<double>.Build.DiagonalIdentity(3, 3));
+
+            var tall = Matrix<double>.Build.Random(8, 3, dist);
+            Assert.IsTrue((tall*Matrix<double>.Build.DiagonalIdentity(3).Multiply(2d)).Equals(tall.Multiply(2d)));
+            Assert.IsTrue((tall*Matrix<double>.Build.Diagonal(3, 5, 2d)).Equals(tall.Multiply(2d).Append(Matrix<double>.Build.Dense(8, 2))));
+            Assert.IsTrue((tall*Matrix<double>.Build.Diagonal(3, 2, 2d)).Equals(tall.Multiply(2d).SubMatrix(0, 8, 0, 2)));
+
+            var wide = Matrix<double>.Build.Random(3, 8, dist);
+            Assert.IsTrue((wide*Matrix<double>.Build.DiagonalIdentity(8).Multiply(2d)).Equals(wide.Multiply(2d)));
+            Assert.IsTrue((wide*Matrix<double>.Build.Diagonal(8, 10, 2d)).Equals(wide.Multiply(2d).Append(Matrix<double>.Build.Dense(3, 2))));
+            Assert.IsTrue((wide*Matrix<double>.Build.Diagonal(8, 2, 2d)).Equals(wide.Multiply(2d).SubMatrix(0, 3, 0, 2)));
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Single/DiagonalMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Single/DiagonalMatrixTests.cs
index 26ae0809..8773fce1 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/DiagonalMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/DiagonalMatrixTests.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using MathNet.Numerics.Distributions;
 using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.LinearAlgebra.Single;
+using MathNet.Numerics.Random;
 using NUnit.Framework;
 
 namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
@@ -377,5 +379,23 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
             var matrix = TestMatrices["Square3x3"];
             Assert.IsTrue(matrix.IsSymmetric);
         }
+
+        [Test]
+        public void DenseDiagonalMatrixMultiplication()
+        {
+            var dist = new ContinuousUniform(-1.0, 1.0, new MersenneTwister());
+
+            Assert.IsInstanceOf<DiagonalMatrix>(Matrix<float>.Build.DiagonalIdentity(3, 3));
+
+            var tall = Matrix<float>.Build.Random(8, 3, dist);
+            Assert.IsTrue((tall*Matrix<float>.Build.DiagonalIdentity(3).Multiply(2f)).Equals(tall.Multiply(2f)));
+            Assert.IsTrue((tall*Matrix<float>.Build.Diagonal(3, 5, 2f)).Equals(tall.Multiply(2f).Append(Matrix<float>.Build.Dense(8, 2))));
+            Assert.IsTrue((tall*Matrix<float>.Build.Diagonal(3, 2, 2f)).Equals(tall.Multiply(2f).SubMatrix(0, 8, 0, 2)));
+
+            var wide = Matrix<float>.Build.Random(3, 8, dist);
+            Assert.IsTrue((wide*Matrix<float>.Build.DiagonalIdentity(8).Multiply(2f)).Equals(wide.Multiply(2f)));
+            Assert.IsTrue((wide*Matrix<float>.Build.Diagonal(8, 10, 2f)).Equals(wide.Multiply(2f).Append(Matrix<float>.Build.Dense(3, 2))));
+            Assert.IsTrue((wide*Matrix<float>.Build.Diagonal(8, 2, 2f)).Equals(wide.Multiply(2f).SubMatrix(0, 3, 0, 2)));
+        }
     }
 }