diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index 19cdf881..04aba072 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -181,5 +181,69 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             Data[column * RowCount + row] = value;
         }
+
+        #region Elementary operations
+        /// <summary>
+        /// Adds another matrix to this matrix. The result will be written into this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to add to this matrix.</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>
+        public override void Add(Matrix other)
+        {
+            if (other == null)
+            {
+                throw new ArgumentNullException("other");
+            }
+
+            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixSameDimensions);
+            }
+
+            DenseMatrix m = other as DenseMatrix;
+            if (m == null)
+            {
+                Parallel.For(0, RowCount, i =>
+                    Parallel.For(0, ColumnCount, j =>
+                        At(i, j, At(i, j) + other.At(i, j))));
+            }
+            else
+            {
+                Control.LinearAlgebraProvider.AddArrays(Data, m.Data, Data);
+            }
+        }
+
+        /// <summary>
+        /// Subtracts another matrix from this matrix. The result will be written into this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to subtract.</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>
+        public override void Subtract(Matrix other)
+        {
+            if (other == null)
+            {
+                throw new ArgumentNullException("other");
+            }
+
+            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixSameDimensions);
+            }
+
+            DenseMatrix m = other as DenseMatrix;
+            if (m == null)
+            {
+                Parallel.For(0, RowCount, i =>
+                    Parallel.For(0, ColumnCount, j =>
+                        At(i, j, At(i, j) - other.At(i, j))));
+            }
+            else
+            {
+                Control.LinearAlgebraProvider.SubtractArrays(Data, m.Data, Data);
+            }
+        }
+        #endregion
     }
 }
diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs
index c5226e2b..e4bcd869 100644
--- a/src/Numerics/LinearAlgebra/Double/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs
@@ -30,6 +30,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 {
     using System;
     using System.Text;
+    using Threading;
 
     using Properties;
 
@@ -411,5 +412,117 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         #endregion
+        
+        #region Elementary operations
+        /// <summary>
+        /// Adds another matrix to this matrix. The result will be written into this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to add to this matrix.</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>
+        public virtual void Add(Matrix other)
+        {
+            if (other == null)
+            {
+                throw new ArgumentNullException("other");
+            }
+
+            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixSameDimensions);
+            }
+
+            Parallel.For(0, RowCount, i => 
+                Parallel.For(0, ColumnCount, j =>
+                    At(i, j, At(i,j) + other.At(i,j))));
+        }
+
+        /// <summary>
+        /// Adds two matrices together and returns the results.
+        /// </summary>
+        /// <remarks>This operator will allocate new memory for the result. It will
+        /// choose the representation of either <paramref name="leftSide"/> or <paramref name="rightSide"/> depending on which
+        /// is denser.</remarks>
+        /// <param name="leftSide">The left matrix to add.</param>
+        /// <param name="rightSide">The right matrix to add.</param>
+        /// <returns>The result of the addition.</returns>
+        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="leftSide"/> and <paramref name="rightSide"/> don't have the same dimensions.</exception>
+        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
+        public static Matrix operator +(Matrix leftSide, Matrix rightSide)
+        {
+            if (rightSide == null)
+            {
+                throw new ArgumentNullException("rightSide");
+            }
+            if (leftSide == null)
+            {
+                throw new ArgumentNullException("leftSide");
+            }
+
+            if (leftSide.RowCount != rightSide.RowCount || leftSide.ColumnCount != rightSide.ColumnCount)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixSameDimensions);
+            }
+
+            Matrix ret = leftSide.Clone();
+            ret.Add(rightSide);
+            return ret;
+        }
+
+        /// <summary>
+        /// Subtracts another matrix from this matrix. The result will be written into this matrix.
+        /// </summary>
+        /// <param name="other">The matrix to subtract.</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>
+        public virtual void Subtract(Matrix other)
+        {
+            if (other == null)
+            {
+                throw new ArgumentNullException("other");
+            }
+
+            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixSameDimensions);
+            }
+
+            Parallel.For(0, RowCount, i =>
+                Parallel.For(0, ColumnCount, j =>
+                    At(i, j, At(i, j) - other.At(i, j))));
+        }
+
+        /// <summary>
+        /// Subtracts two matrices together and returns the results.
+        /// </summary>
+        /// <remarks>This operator will allocate new memory for the result. It will
+        /// choose the representation of either <paramref name="leftSide"/> or <paramref name="rightSide"/> depending on which
+        /// is denser.</remarks>
+        /// <param name="leftSide">The left matrix to subtract.</param>
+        /// <param name="rightSide">The right matrix to subtract.</param>
+        /// <returns>The result of the addition.</returns>
+        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="leftSide"/> and <paramref name="rightSide"/> don't have the same dimensions.</exception>
+        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
+        public static Matrix operator -(Matrix leftSide, Matrix rightSide)
+        {
+            if (rightSide == null)
+            {
+                throw new ArgumentNullException("rightSide");
+            }
+            if (leftSide == null)
+            {
+                throw new ArgumentNullException("leftSide");
+            }
+
+            if (leftSide.RowCount != rightSide.RowCount || leftSide.ColumnCount != rightSide.ColumnCount)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixSameDimensions);
+            }
+
+            Matrix ret = leftSide.Clone();
+            ret.Subtract(rightSide);
+            return ret;
+        }
+        #endregion
     }
 }
\ No newline at end of file
diff --git a/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs
index 3bd95a35..ff88a696 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs
@@ -7,10 +7,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
 
     public class DenseMatrixTests : MatrixTests
     {
-        protected double[] SquareMatrix = new double[] { -1.1, 0.0, 1.0, -4.4, -2.2, 1.1, 2.1, 5.5, -3.3, 2.2, 6.2, 6.6, -4.4, 3.3, 4.3, -7.7 };
-        protected int SquareMatrixRows = 4;
-        protected int SquareMatrixColumns = 4;
-
         protected override Matrix CreateMatrix(int rows, int columns)
         {
             return new DenseMatrix(rows, columns);
diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs
index 41ed50a7..0b81e099 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs
@@ -15,12 +15,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         protected abstract Matrix CreateMatrix(double[,] data);
 
         [SetUp]
-        public void SetupDistributions()
+        public void SetupMatrices()
         {
             testData2D = new Dictionary<string, double[,]>();
             testData2D.Add("Singular3x3", new double[,] { { 1, 1, 2 }, { 1, 1, 2 }, { 1, 1, 2 } });
             testData2D.Add("Square3x3", new double[,] { { -1.1, -2.2, -3.3 }, { 0, 1.1, 2.2 }, { -4.4, 5.5, 6.6 } });
             testData2D.Add("Square4x4", new double[,] { { -1.1, -2.2, -3.3, -4.4 }, { 0, 1.1, 2.2, 3.3 }, { 1.0, 2.1, 6.2, 4.3 }, { -4.4, 5.5, 6.6, -7.7 } });
+            testData2D.Add("Singular4x4", new double[,] { { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 } });
             testData2D.Add("Tall3x2", new double[,] { { -1.1, -2.2 }, { 0, 1.1 }, { -4.4, 5.5 } });
             testData2D.Add("Wide2x3", new double[,] { { -1.1, -2.2, -3.3 }, { 0, 1.1, 2.2 } });
 
@@ -160,5 +161,207 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         public void MatrixGetHashCode()
         {
         }
+
+        #region Elementary operations
+        [Test]
+        [Row("Singular3x3", "Square3x3")]
+        [Row("Singular4x4", "Square4x4")]
+        public void AddMatrix(string mtxA, string mtxB)
+        {
+            var A = testMatrices[mtxA];
+            var B = testMatrices[mtxB];
+
+            Matrix matrix = A.Clone();
+            matrix.Add(B);
+            for (int i = 0; i < matrix.RowCount; i++)
+            {
+                for (int j = 0; j < matrix.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrix[i, j], A[i, j] + B[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void AddMatrixThrowsExceptionWhenArgumentIsNull()
+        {
+            Matrix matrix = testMatrices["Singular4x4"];
+            Matrix other = null;
+            matrix.Add(other);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void AddMatrixThrowsExceptionArgumentHasTooFewColumns()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Tall3x2"];
+            matrix.Add(other);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void AddMatrixThrowsExceptionArgumentHasTooFewRows()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Wide2x3"];
+            matrix.Add(other);
+        }
+
+        [Test]
+        [Row("Singular3x3", "Square3x3")]
+        [Row("Singular4x4", "Square4x4")]
+        public void AddOperator(string mtxA, string mtxB)
+        {
+            var A = testMatrices[mtxA];
+            var B = testMatrices[mtxB];
+            
+            Matrix result = A + B;
+            for (int i = 0; i < A.RowCount; i++)
+            {
+                for (int j = 0; j < A.ColumnCount; j++)
+                {
+                    Assert.AreEqual(result[i,j], A[i, j] + B[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void AddOperatorThrowsExceptionWhenLeftsideIsNull()
+        {
+            Matrix matrix = null;
+            Matrix other = testMatrices["Singular3x3"];
+            Matrix result = matrix + other;
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void AddOperatorThrowsExceptionWhenRightsideIsNull()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = null;
+            Matrix result = matrix + other;
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void AddOperatorThrowsExceptionWhenRightsideHasTooFewColumns()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Tall3x2"];
+            Matrix result = matrix + other;
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void AddOperatorThrowsExceptionWhenRightsideHasTooFewRows()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Wide2x3"];
+            Matrix result = matrix + other;
+        }
+
+        [Test]
+        [Row("Singular3x3", "Square3x3")]
+        [Row("Singular4x4", "Square4x4")]
+        public void SubtractMatrix(string mtxA, string mtxB)
+        {
+            var A = testMatrices[mtxA];
+            var B = testMatrices[mtxB];
+
+            Matrix matrix = A.Clone();
+            matrix.Subtract(B);
+            for (int i = 0; i < matrix.RowCount; i++)
+            {
+                for (int j = 0; j < matrix.ColumnCount; j++)
+                {
+                    Assert.AreEqual(matrix[i, j], A[i, j] - B[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void SubtractMatrixThrowsExceptionWhenRightSideIsNull()
+        {
+            Matrix matrix = testMatrices["Singular4x4"];
+            Matrix other = null;
+            matrix.Subtract(other);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void SubtractMatrixThrowsExceptionWhenRightSideHasTooFewColumns()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Tall3x2"];
+            matrix.Subtract(other);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void SubtractMatrixThrowsExceptionWhenRightSideHasTooFewRows()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Wide2x3"];
+            matrix.Subtract(other);
+        }
+
+        [Test]
+        [Row("Singular3x3", "Square3x3")]
+        [Row("Singular4x4", "Square4x4")]
+        public void SubtractOperator(string mtxA, string mtxB)
+        {
+            var A = testMatrices[mtxA];
+            var B = testMatrices[mtxB];
+
+            Matrix result = A - B;
+            for (int i = 0; i < A.RowCount; i++)
+            {
+                for (int j = 0; j < A.ColumnCount; j++)
+                {
+                    Assert.AreEqual(result[i, j], A[i, j] - B[i, j]);
+                }
+            }
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void SubtractOperatorThrowsExceptionWhenLeftsideIsNull()
+        {
+            Matrix matrix = null;
+            Matrix other = testMatrices["Singular3x3"];
+            Matrix result = matrix - other;
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void SubtractOperatorThrowsExceptionWhenRightsideIsNull()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = null;
+            Matrix result = matrix - other;
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void SubtractOperatorThrowsExceptionWhenRightsideHasTooFewColumns()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Tall3x2"];
+            Matrix result = matrix - other;
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void SubtractOperatorThrowsExceptionWhenRightsideHasTooFewRows()
+        {
+            Matrix matrix = testMatrices["Singular3x3"];
+            Matrix other = testMatrices["Wide2x3"];
+            Matrix result = matrix - other;
+        }
+        #endregion
     }
 }
\ No newline at end of file