diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index c8bf954b..f240e1d5 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
     using Algorithms.LinearAlgebra;
     using Generic;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// A Matrix class with dense storage. The underlying storage is a one dimensional array in column-major order.
@@ -240,50 +239,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, RowCount, ColumnCount, Data);
         }
 
-        #region Elementary operations
-
-        /// <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(Data, denseOther.Data, denseResult.Data);
-            }
-        }
-
-        /// <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(Data, denseOther.Data, denseResult.Data);
-            }
-        }
-    
-        #endregion
-
         #region Static constructors for special matrices.
 
         /// <summary>
@@ -307,25 +262,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
         #endregion
 
-        /// <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.Data[(i * ColumnCount) + j] = Data[index + i].Conjugate();
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>
         /// Multiplies each element of the matrix by a scalar and places results into the result matrix.
         /// </summary>
@@ -526,6 +462,65 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
+        /// <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(Data, denseOther.Data, denseResult.Data);
+            }
+        }
+
+        /// <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(Data, denseOther.Data, denseResult.Data);
+            }
+        }
+
+        /// <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.Data[(i * ColumnCount) + j] = Data[index + i].Conjugate();
+                }
+            }
+
+            return ret;
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
index 657daa2e..a3db6c85 100644
--- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
@@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
     using Distributions;
     using Generic;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// <c>Complex</c> version of the <see cref="Matrix{T}"/> class.
@@ -188,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
-         /// <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>
@@ -327,6 +326,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(Complex divisor, Matrix<Complex> result)
+        {
+            throw new NotImplementedException();
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 1332cc55..bd948a02 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
     using Generic;
     using Numerics;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// A Matrix class with dense storage. The underlying storage is a one dimensional array in column-major order.
@@ -240,50 +239,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, RowCount, ColumnCount, Data);
         }
 
-        #region Elementary operations
-
-        /// <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(Data, denseOther.Data, denseResult.Data);
-            }
-        }
-
-        /// <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(Data, denseOther.Data, denseResult.Data);
-            }
-        }
-    
-        #endregion
-
         #region Static constructors for special matrices.
 
         /// <summary>
@@ -307,25 +262,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
         #endregion
 
-        /// <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.Data[(i * ColumnCount) + j] = Data[index + i].Conjugate();
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>
         /// Multiplies each element of the matrix by a scalar and places results into the result matrix.
         /// </summary>
@@ -526,6 +462,65 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
+        /// <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(Data, denseOther.Data, denseResult.Data);
+            }
+        }
+
+        /// <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(Data, denseOther.Data, denseResult.Data);
+            }
+        }
+
+        /// <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.Data[(i * ColumnCount) + j] = Data[index + i].Conjugate();
+                }
+            }
+
+            return ret;
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
index 833f12b3..0fbe7008 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
@@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
     using Generic;
     using Numerics;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// <c>Complex32</c> version of the <see cref="Matrix{T}"/> class.
@@ -327,6 +326,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(Complex32 divisor, Matrix<Complex32> result)
+        {
+            throw new NotImplementedException();
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index 73332c9c..800cec08 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -239,7 +239,28 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, RowCount, ColumnCount, Data);
         }
 
-        #region Elementary operations
+        #region Static constructors for special matrices.
+
+        /// <summary>
+        /// Initializes a square <see cref="DenseMatrix"/> with all zero's except for ones on the diagonal.
+        /// </summary>
+        /// <param name="order">the size of the square matrix.</param>
+        /// <returns>A dense identity matrix.</returns>
+        /// <exception cref="ArgumentException">
+        /// If <paramref name="order"/> is less than one.
+        /// </exception>
+        public static DenseMatrix Identity(int order)
+        {
+            var m = new DenseMatrix(order);
+            for (var i = 0; i < order; i++)
+            {
+                m.Data[(i * order) + i] = 1.0;
+            }
+
+            return m;
+        }
+
+        #endregion
 
         /// <summary>
         /// Adds another matrix to this matrix.
@@ -281,39 +302,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
     
-        #endregion
-
-        #region Static constructors for special matrices.
-
-        /// <summary>
-        /// Initializes a square <see cref="DenseMatrix"/> with all zero's except for ones on the diagonal.
-        /// </summary>
-        /// <param name="order">the size of the square matrix.</param>
-        /// <returns>A dense identity matrix.</returns>
-        /// <exception cref="ArgumentException">
-        /// If <paramref name="order"/> is less than one.
-        /// </exception>
-        public static DenseMatrix Identity(int order)
-        {
-            var m = new DenseMatrix(order);
-            for (var i = 0; i < order; i++)
-            {
-                m.Data[(i * order) + i] = 1.0;
-            }
-
-            return m;
-        }
-
-        #endregion
-
-        /// <summary>
-        /// Returns the conjugate transpose of this matrix.
-        /// </summary>        
-        /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<double> ConjugateTranspose()
-        {
-            return Transpose();
-        }
 
         /// <summary>
         /// Multiplies each element of the matrix by a scalar and places results into the result matrix.
@@ -515,6 +503,42 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(double divisor, Matrix<double> result)
+        {
+            var denseResult = result as DenseMatrix;
+
+            if (denseResult == null)
+            {
+               base.DoModulus(divisor, result);
+            }
+            else
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
+                }
+
+                CommonParallel.For(
+                    0, 
+                    Data.Length,
+                    index => denseResult.Data[index] %= divisor);
+            }
+        }
+
+        /// <summary>
+        /// Returns the conjugate transpose of this matrix.
+        /// </summary>        
+        /// <returns>The conjugate transpose of this matrix.</returns>
+        public override Matrix<double> ConjugateTranspose()
+        {
+            return Transpose();
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
index 18a51154..5550b5aa 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
@@ -632,7 +632,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return (DenseVector)leftSide.Modulus(rightSide);
         }
 
-
         /// <summary>
         /// Returns the index of the absolute minimum element.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 47a86955..ffcdb900 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -1547,6 +1547,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(double divisor, Matrix<double> result)
+        {
+            var denseResult = result as DiagonalMatrix;
+
+            if (denseResult == null)
+            {
+                base.DoModulus(divisor, result);
+            }
+            else
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
+                }
+
+                for (var index = 0; index < Data.Length; index++)
+                {
+                    denseResult.Data[index] %= divisor;
+                }
+            }
+        }
+
         #region Static constructors for special matrices.
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs
index 65a0aace..0e4e546e 100644
--- a/src/Numerics/LinearAlgebra/Double/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs
@@ -30,7 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
     using Distributions;
     using Generic;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// <c>double</c> version of the <see cref="Matrix{T}"/> class.
@@ -317,6 +316,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(double divisor, Matrix<double> result)
+        {
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    result.At(row, column, At(row, column) % divisor);
+                }
+            }
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index a7c7a502..8d197d7c 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
@@ -1483,6 +1483,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(double divisor, Matrix<double> result)
+        {
+            var denseResult = result as SparseMatrix;
+
+            if (denseResult == null)
+            {
+                base.DoModulus(divisor, result);
+            }
+            else
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
+                }
+
+                for (var index = 0; index < denseResult._nonZeroValues.Length; index++)
+                {
+                    denseResult._nonZeroValues[index] %= divisor;
+                }
+            }
+        }
+
         /// <summary>
         /// Iterates throw each element in the matrix (row-wise).
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs
index 51ab1f3b..c0ab7de3 100644
--- a/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs
+++ b/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs
@@ -30,7 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
     using Distributions;
     using Factorization;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// Defines the base class for <c>Matrix</c> classes.
@@ -218,7 +217,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
                 return;
             }
 
-            //CopyTo(result);
             DoMultiply(scalar, result);
         }
 
@@ -933,6 +931,62 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <param name="result">The matrix to store the result of the pointwise division.</param>
         protected abstract void DoPointwiseDivide(Matrix<T> other, Matrix<T> result);
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <returns>A matrix containing the results.</returns>
+        public Matrix<T> Modulus(T divisor)
+        {
+            var result = CreateMatrix(RowCount, ColumnCount);
+            DoModulus(divisor, result);
+            return result;
+        }
+
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        public void Modulus(T divisor, Matrix<T> result)
+        {
+            if (result == null)
+            {
+                throw new ArgumentNullException("result");
+            }
+
+            if (ColumnCount != result.ColumnCount || RowCount != result.RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result");
+            }
+
+            DoModulus(divisor, result);
+        }
+
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected abstract void DoModulus(T divisor, Matrix<T> result);
+
+        /// <summary>
+        /// Multiplies a <strong>Matrix</strong> by a constant and returns the result.
+        /// </summary>
+        /// <param name="leftSide">The matrix to multiply.</param>
+        /// <param name="rightSide">The constant to multiply the matrix by.</param>
+        /// <returns>The result of the multiplication.</returns>
+        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> is <see langword="null" />.</exception>
+        public static Matrix<T> operator %(Matrix<T> leftSide, T rightSide)
+        {
+            if (leftSide == null)
+            {
+                throw new ArgumentNullException("leftSide");
+            }
+
+            return leftSide.Modulus(rightSide);
+        }
+
         /// <summary>
         /// Generates a matrix with random elements.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 0b4d0027..3a6ca20b 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -239,7 +239,28 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, RowCount, ColumnCount, Data);
         }
 
-        #region Elementary operations
+        #region Static constructors for special matrices.
+
+        /// <summary>
+        /// Initializes a square <see cref="DenseMatrix"/> with all zero's except for ones on the diagonal.
+        /// </summary>
+        /// <param name="order">the size of the square matrix.</param>
+        /// <returns>A dense identity matrix.</returns>
+        /// <exception cref="ArgumentException">
+        /// If <paramref name="order"/> is less than one.
+        /// </exception>
+        public static DenseMatrix Identity(int order)
+        {
+            var m = new DenseMatrix(order);
+            for (var i = 0; i < order; i++)
+            {
+                m.Data[(i * order) + i] = 1.0f;
+            }
+
+            return m;
+        }
+
+        #endregion
 
         /// <summary>
         /// Adds another matrix to this matrix.
@@ -280,40 +301,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 Control.LinearAlgebraProvider.SubtractArrays(Data, denseOther.Data, denseResult.Data);
             }
         }
-    
-        #endregion
-
-        #region Static constructors for special matrices.
-
-        /// <summary>
-        /// Initializes a square <see cref="DenseMatrix"/> with all zero's except for ones on the diagonal.
-        /// </summary>
-        /// <param name="order">the size of the square matrix.</param>
-        /// <returns>A dense identity matrix.</returns>
-        /// <exception cref="ArgumentException">
-        /// If <paramref name="order"/> is less than one.
-        /// </exception>
-        public static DenseMatrix Identity(int order)
-        {
-            var m = new DenseMatrix(order);
-            for (var i = 0; i < order; i++)
-            {
-                m.Data[(i * order) + i] = 1.0f;
-            }
-
-            return m;
-        }
-
-        #endregion
-
-        /// <summary>
-        /// Returns the conjugate transpose of this matrix.
-        /// </summary>        
-        /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<float> ConjugateTranspose()
-        {
-            return Transpose();
-        }
 
         /// <summary>
         /// Multiplies each element of the matrix by a scalar and places results into the result matrix.
@@ -515,6 +502,42 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(float divisor, Matrix<float> result)
+        {
+            var denseResult = result as DenseMatrix;
+
+            if (denseResult == null)
+            {
+                base.DoModulus(divisor, result);
+            }
+            else
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
+                }
+
+                CommonParallel.For(
+                    0,
+                    Data.Length,
+                    index => denseResult.Data[index] %= divisor);
+            }
+        }
+
+        /// <summary>
+        /// Returns the conjugate transpose of this matrix.
+        /// </summary>        
+        /// <returns>The conjugate transpose of this matrix.</returns>
+        public override Matrix<float> ConjugateTranspose()
+        {
+            return Transpose();
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index 72c93ae6..5bac5ed6 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -1547,6 +1547,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(float divisor, Matrix<float> result)
+        {
+            var denseResult = result as DiagonalMatrix;
+
+            if (denseResult == null)
+            {
+                base.DoModulus(divisor, result);
+            }
+            else
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
+                }
+
+                for (var index = 0; index < Data.Length; index++)
+                {
+                    denseResult.Data[index] %= divisor;
+                }
+            }
+        }
+
         #region Static constructors for special matrices.
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs
index 998046c5..3214c33e 100644
--- a/src/Numerics/LinearAlgebra/Single/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs
@@ -30,7 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
     using Distributions;
     using Generic;
     using Properties;
-    using Threading;
 
     /// <summary>
     /// <c>float</c> version of the <see cref="Matrix{T}"/> class.
@@ -270,6 +269,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(float divisor, Matrix<float> result)
+        {
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    result.At(row, column, At(row, column) % divisor);
+                }
+            }
+        }
+
         /// <summary>
         /// Negate each element of this matrix and place the results into the result matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index 5eadb4e3..9fa5c8d0 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
@@ -1466,6 +1466,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        /// <summary>
+        /// Computes the modulus for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoModulus(float divisor, Matrix<float> result)
+        {
+            var denseResult = result as SparseMatrix;
+
+            if (denseResult == null)
+            {
+                base.DoModulus(divisor, result);
+            }
+            else
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
+                }
+
+                for (var index = 0; index < denseResult._nonZeroValues.Length; index++)
+                {
+                    denseResult._nonZeroValues[index] %= divisor;
+                }
+            }
+        }
+
         /// <summary>
         /// Iterates throw each element in the matrix (row-wise).
         /// </summary>
diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
index 75f35ba4..e51816df 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
@@ -1001,5 +1001,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             var matrix = TestMatrices["Wide2x3"];
             Assert.Throws<ArgumentException>(() => matrix.Trace());
         }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector.
+        /// </summary>
+        [Test]
+        public void CanComputeModulus()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = matrix.Modulus(3.2);
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2, mod[row, column], 14);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector using a result vector.
+        /// </summary>
+        [Test]
+        public void CanComputeModulusUsingResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
+            matrix.Modulus(3.2, mod);
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2, mod[row, column], 14);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector using a result vector.
+        /// </summary>
+        [Test]
+        public void CanComputeModulusUsingSameResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"].Clone();
+            matrix.Modulus(3.2, matrix);
+            var data = TestMatrices["Square3x3"];
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(data[row, column] % 3.2, matrix[row, column], 14);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector using the operator %.
+        /// </summary>
+        [Test]
+        public void CanComputeModulusUsingOperator()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = matrix % 3.2;
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2, mod[row, column], 14);
+                }
+            }
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs
index 80441647..a07ab17b 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs
@@ -1081,39 +1081,41 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         public void CanComputeModulus()
         {
             var vector = CreateVector(Data);
-            var mod = vector.Modulus(3.2f);
+            var mod = vector.Modulus(3.2);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqual(Data[index] % 3.2, mod[index], 7);
+                AssertHelpers.AlmostEqual(Data[index] % 3.2, mod[index], 14);
             }
         }
 
         /// <summary>
         /// Can compute the modules of each element of vector using a result vector.
         /// </summary>
+        [Test]
         public void CanComputeModulusUsingResultVector()
         {
             var vector = CreateVector(Data);
             var mod = CreateVector(vector.Count);
-            vector.Modulus(3.2f, mod);
+            vector.Modulus(3.2, mod);
 
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqual(Data[index] % 3.2, mod[index], 7);
+                AssertHelpers.AlmostEqual(Data[index] % 3.2, mod[index], 14);
             }
         }
 
         /// <summary>
         /// Can compute the modules of each element of vector using a result vector.
         /// </summary>
+        [Test]
         public void CanComputeModulusUsingSameResultVector()
         {
             var vector = CreateVector(Data);
-            vector.Modulus(3.2f, vector);
+            vector.Modulus(3.2, vector);
 
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqual(Data[index] % 3.2, vector[index], 7);
+                AssertHelpers.AlmostEqual(Data[index] % 3.2, vector[index], 14);
             }
         }
 
@@ -1124,10 +1126,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         public void CanComputeModulusUsingOperator()
         {
             var vector = CreateVector(Data);
-            var mod = vector % 4.5f;
+            var mod = vector % 4.5;
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqual(Data[index] % 4.5, mod[index], 7);
+                AssertHelpers.AlmostEqual(Data[index] % 4.5, mod[index], 14);
             }
         }
     }
diff --git a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs
index 02e3ec43..fbb99b45 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs
@@ -1001,5 +1001,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
             var matrix = TestMatrices["Wide2x3"];
             Assert.Throws<ArgumentException>(() => matrix.Trace());
         }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector.
+        /// </summary>
+        [Test]
+        public void CanComputeModulus()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = matrix.Modulus(3.2f);
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2f, mod[row, column], 7);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector using a result vector.
+        /// </summary>
+        [Test]
+        public void CanComputeModulusUsingResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
+            matrix.Modulus(3.2f, mod);
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2f, mod[row, column], 7);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector using a result vector.
+        /// </summary>
+        [Test]
+        public void CanComputeModulusUsingSameResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"].Clone();
+            matrix.Modulus(3.2f, matrix);
+            var data = TestMatrices["Square3x3"];
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(data[row, column] % 3.2f, matrix[row, column], 7);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can compute the modules of each element of vector using the operator %.
+        /// </summary>
+        [Test]
+        public void CanComputeModulusUsingOperator()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = matrix % 3.2f;
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2f, mod[row, column], 7);
+                }
+            }
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs
index 4d02f587..2528c7d7 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs
@@ -1091,6 +1091,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
         /// <summary>
         /// Can compute the modules of each element of vector using a result vector.
         /// </summary>
+        [Test]
         public void CanComputeModulusUsingResultVector()
         {
             var vector = CreateVector(Data);
@@ -1106,6 +1107,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
         /// <summary>
         /// Can compute the modules of each element of vector using a result vector.
         /// </summary>
+        [Test]
         public void CanComputeModulusUsingSameResultVector()
         {
             var vector = CreateVector(Data);