diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.cs
index 876f5332..4f65922a 100644
--- a/src/Numerics/LinearAlgebra/Generic/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Generic/Matrix.cs
@@ -35,8 +35,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
     using Numerics;
     using Properties;
     using Storage;
-    using Threading;
-
 
     /// <summary>
     /// Defines the base class for <c>Matrix</c> classes.
@@ -290,6 +288,19 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             return result;
         }
 
+#if !PORTABLE
+        /// <summary>
+        /// Creates a new object that is a copy of the current instance.
+        /// </summary>
+        /// <returns>
+        /// A new object that is a copy of this instance.
+        /// </returns>
+        object ICloneable.Clone()
+        {
+            return Clone();
+        }
+#endif
+
         /// <summary>
         /// Copies the elements of this matrix to the given matrix.
         /// </summary>
@@ -339,17 +350,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </remarks>
         public abstract Vector<T> CreateVector(int size, bool fullyMutable = false);
 
-        /// <summary>
-        /// Returns a <see cref="System.String"/> that represents this instance.
-        /// </summary>
-        /// <returns>
-        /// A <see cref="System.String"/> that represents this instance.
-        /// </returns>
-        public override string ToString()
-        {
-            return ToString(null);
-        }
-
         /// <summary>
         /// Copies a row into an Vector.
         /// </summary>
@@ -400,7 +400,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// or greater than or equal to the number of rows.</item>
         /// <item><paramref name="columnIndex"/> is negative,
         /// or greater than or equal to the number of columns.</item>
-        /// <item><c>(columnIndex + length) &gt;= Columns.</c></item></list></exception>        
+        /// <item><c>(columnIndex + length) &gt;= Columns.</c></item></list></exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="length"/> is not positive.</exception>
         public Vector<T> Row(int rowIndex, int columnIndex, int length)
         {
@@ -418,10 +418,10 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <param name="result">The Vector to copy the column into.</param>
         /// <exception cref="ArgumentNullException">If the result Vector is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> is negative,
-        /// or greater than or equal to the number of columns.</exception>        
+        /// or greater than or equal to the number of columns.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is negative,
-        /// or greater than or equal to the number of rows.</exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> + <paramref name="length"/>  
+        /// or greater than or equal to the number of rows.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> + <paramref name="length"/>
         /// is greater than or equal to the number of rows.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="length"/> is not positive.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <strong>result.Count &lt; length</strong>.</exception>
@@ -486,7 +486,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <item><paramref name="rowIndex"/> is negative,
         /// or greater than or equal to the number of rows.</item>
         /// <item><c>(rowIndex + length) &gt;= Rows.</c></item></list>
-        /// </exception>        
+        /// </exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="length"/> is not positive.</exception>
         public Vector<T> Column(int columnIndex, int rowIndex, int length)
         {
@@ -504,10 +504,10 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <param name="result">The Vector to copy the column into.</param>
         /// <exception cref="ArgumentNullException">If the result Vector is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is negative,
-        /// or greater than or equal to the number of columns.</exception>        
+        /// or greater than or equal to the number of columns.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> is negative,
-        /// or greater than or equal to the number of rows.</exception>        
-        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> + <paramref name="length"/>  
+        /// or greater than or equal to the number of rows.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> + <paramref name="length"/>
         /// is greater than or equal to the number of rows.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="length"/> is not positive.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <strong>result.Count &lt; length</strong>.</exception>
@@ -521,10 +521,29 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             Storage.CopySubColumnTo(result.Storage, columnIndex, rowIndex, 0, length);
         }
 
+        /// <summary>
+        /// Returns a new matrix containing the upper triangle of this matrix.
+        /// </summary>
+        /// <returns>The upper triangle of this matrix.</returns>
+        public virtual Matrix<T> UpperTriangle()
+        {
+            var ret = CreateMatrix(RowCount, ColumnCount);
+
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = row; column < ColumnCount; column++)
+                {
+                    ret.At(row, column, At(row, column));
+                }
+            }
+
+            return ret;
+        }
+
         /// <summary>
         /// Returns a new matrix containing the lower triangle of this matrix.
         /// </summary>
-        /// <returns>The lower triangle of this matrix.</returns>        
+        /// <returns>The lower triangle of this matrix.</returns>
         public virtual Matrix<T> LowerTriangle()
         {
             var ret = CreateMatrix(RowCount, ColumnCount);
@@ -567,25 +586,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             }
         }
 
-        /// <summary>
-        /// Returns a new matrix containing the upper triangle of this matrix.
-        /// </summary>
-        /// <returns>The upper triangle of this matrix.</returns>   
-        public virtual Matrix<T> UpperTriangle()
-        {
-            var ret = CreateMatrix(RowCount, ColumnCount);
-
-            for (var row = 0; row < RowCount; row++)
-            {
-                for (var column = row; column < ColumnCount; column++)
-                {
-                    ret.At(row, column, At(row, column));
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>
         /// Puts the upper triangle of this matrix into the result matrix.
         /// </summary>
@@ -623,10 +623,10 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <returns>The requested sub-matrix.</returns>
         /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
         /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
+        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number
         /// of columns.</item>
         /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>        
+        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
         /// is not positive.</exception>
         public virtual Matrix<T> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
@@ -636,105 +636,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             return target;
         }
 
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates over the matrix columns.
-        /// </summary>
-        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over the matrix columns</returns>
-        /// <seealso cref="IEnumerator{T}"/>
-        public virtual IEnumerable<Tuple<int, Vector<T>>> ColumnEnumerator()
-        {
-            for (var i = 0; i < ColumnCount; i++)
-            {
-                yield return new Tuple<int, Vector<T>>(i, Column(i));
-            }
-        }
-
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates the requested matrix columns.
-        /// </summary>
-        /// <param name="index">The column to start enumerating over.</param>
-        /// <param name="length">The number of columns to enumerating over.</param>
-        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over requested matrix columns.</returns>
-        /// <seealso cref="IEnumerator{T}"/>
-        /// <exception cref="ArgumentOutOfRangeException">If:
-        /// <list><item><paramref name="index"/> is negative,
-        /// or greater than or equal to the number of columns.</item>        
-        /// <item><c>(index + length) &gt;= Columns.</c></item></list>
-        /// </exception>   
-        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>     
-        public virtual IEnumerable<Tuple<int, Vector<T>>> ColumnEnumerator(int index, int length)
-        {
-            if (index >= ColumnCount || index < 0)
-            {
-                throw new ArgumentOutOfRangeException("index");
-            }
-
-            if (index + length > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("length");
-            }
-
-            if (length < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
-            }
-
-            var maxIndex = index + length;
-            for (var i = index; i < maxIndex; i++)
-            {
-                yield return new Tuple<int, Vector<T>>(i, Column(i));
-            }
-        }
-
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates the requested matrix rows.
-        /// </summary>
-        /// <param name="index">The row to start enumerating over.</param>
-        /// <param name="length">The number of rows to enumerating over.</param>
-        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over requested matrix rows.</returns>
-        /// <seealso cref="IEnumerator{T}"/>
-        /// <exception cref="ArgumentOutOfRangeException">If:
-        /// <list><item><paramref name="index"/> is negative,
-        /// or greater than or equal to the number of rows.</item>        
-        /// <item><c>(index + length) &gt;= Rows.</c></item></list></exception>        
-        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>     
-        public virtual IEnumerable<Tuple<int, Vector<T>>> RowEnumerator(int index, int length)
-        {
-            if (index >= RowCount || index < 0)
-            {
-                throw new ArgumentOutOfRangeException("index");
-            }
-
-            if (index + length > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("length");
-            }
-
-            if (length < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
-            }
-
-            var maxi = index + length;
-            for (var i = index; i < maxi; i++)
-            {
-                yield return new Tuple<int, Vector<T>>(i, Row(i));
-            }
-        }
-
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates over the matrix rows.
-        /// </summary>
-        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over the matrix rows</returns>
-        /// <seealso cref="IEnumerator{T}"/>        
-        public virtual IEnumerable<Tuple<int, Vector<T>>> RowEnumerator()
-        {
-            for (var i = 0; i < RowCount; i++)
-            {
-                yield return new Tuple<int, Vector<T>>(i, Row(i));
-            }
-        }
-
         /// <summary>
         /// Returns the elements of the diagonal in a Vector.
         /// </summary>
@@ -896,7 +797,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="columnIndex">The column to copy the values to.</param>
         /// <param name="column">The array to copy the values from.</param>
-        /// <exception cref="ArgumentNullException">If <paramref name="column"/> is <see langword="null" />.</exception>        
+        /// <exception cref="ArgumentNullException">If <paramref name="column"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is less than zero,
         /// or greater than or equal to the number of columns.</exception>
         /// <exception cref="ArgumentException">If the size of <paramref name="column"/> does not
@@ -931,7 +832,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="columnIndex">The column to copy the values to.</param>
         /// <param name="column">The vector to copy the values from.</param>
-        /// <exception cref="ArgumentNullException">If <paramref name="column"/> is <see langword="null" />.</exception>        
+        /// <exception cref="ArgumentNullException">If <paramref name="column"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is less than zero,
         /// or greater than or equal to the number of columns.</exception>
         /// <exception cref="ArgumentException">If the size of <paramref name="column"/> does not
@@ -1007,7 +908,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="rowIndex">The row to copy the values to.</param>
         /// <param name="row">The vector to copy the values from.</param>
-        /// <exception cref="ArgumentNullException">If <paramref name="row"/> is <see langword="null" />.</exception>            
+        /// <exception cref="ArgumentNullException">If <paramref name="row"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> is less than zero,
         /// or greater than or equal to the number of rows.</exception>
         /// <exception cref="ArgumentException">If the size of <paramref name="row"/> does not
@@ -1040,7 +941,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="rowIndex">The row to copy the values to.</param>
         /// <param name="row">The array to copy the values from.</param>
-        /// <exception cref="ArgumentNullException">If <paramref name="row"/> is <see langword="null" />.</exception>  
+        /// <exception cref="ArgumentNullException">If <paramref name="row"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> is less than zero,
         /// or greater than or equal to the number of rows.</exception>
         /// <exception cref="ArgumentException">If the size of <paramref name="row"/> does not
@@ -1078,10 +979,10 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <param name="subMatrix">The sub-matrix to copy from.</param>
         /// <exception cref="ArgumentOutOfRangeException">If: <list><item><paramref name="rowIndex"/> is
         /// negative, or greater than or equal to the number of rows.</item>
-        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number 
+        /// <item><paramref name="columnIndex"/> is negative, or greater than or equal to the number
         /// of columns.</item>
         /// <item><c>(columnIndex + columnLength) &gt;= Columns</c></item>
-        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception> 
+        /// <item><c>(rowIndex + rowLength) &gt;= Rows</c></item></list></exception>
         /// <exception cref="ArgumentNullException">If <paramref name="subMatrix"/> is <see langword="null" /></exception>
         /// <item>the size of <paramref name="subMatrix"/> is not at least <paramref name="rowCount"/> x <paramref name="columnCount"/>.</item>
         /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowCount"/> or <paramref name="columnCount"/>
@@ -1101,7 +1002,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="source">The vector to copy the values from. The length of the vector should be
         /// Min(Rows, Columns).</param>
-        /// <exception cref="ArgumentNullException">If <paramref name="source"/> is <see langword="null" />.</exception>   
+        /// <exception cref="ArgumentNullException">If <paramref name="source"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If the length of <paramref name="source"/> does not
         /// equal Min(Rows, Columns).</exception>
         /// <remarks>For non-square matrices, the elements of <paramref name="source"/> are copied to
@@ -1131,7 +1032,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <param name="source">The array to copy the values from. The length of the vector should be
         /// Min(Rows, Columns).</param>
-        /// <exception cref="ArgumentNullException">If <paramref name="source"/> is <see langword="null" />.</exception>   
+        /// <exception cref="ArgumentNullException">If <paramref name="source"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If the length of <paramref name="source"/> does not
         /// equal Min(Rows, Columns).</exception>
         /// <remarks>For non-square matrices, the elements of <paramref name="source"/> are copied to
@@ -1157,221 +1058,71 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         }
 
         /// <summary>
-        /// Returns this matrix as a multidimensional array.
+        /// Returns the transpose of this matrix.
         /// </summary>
-        /// <returns>A multidimensional containing the values of this matrix.</returns>        
-        public T[,] ToArray()
+        /// <returns>The transpose of this matrix.</returns>
+        public virtual Matrix<T> Transpose()
         {
-            return Storage.ToArray();
-        }
+            var ret = CreateMatrix(ColumnCount, RowCount);
+            for (var j = 0; j < ColumnCount; j++)
+            {
+                for (var i = 0; i < RowCount; i++)
+                {
+                    ret.At(j, i, At(i, j));
+                }
+            }
 
-        /// <summary>
-        /// Returns the matrix's elements as an array with the data laid out column-wise.
-        /// </summary>
-        /// <example><pre>
-        /// 1, 2, 3
-        /// 4, 5, 6  will be returned as  1, 4, 7, 2, 5, 8, 3, 6, 9
-        /// 7, 8, 9
-        /// </pre></example>
-        /// <returns>An array containing the matrix's elements.</returns>
-        public T[] ToColumnWiseArray()
-        {
-            return Storage.ToColumnMajorArray();
+            return ret;
         }
 
         /// <summary>
-        /// Returns the matrix's elements as an array with the data laid row-wise.
+        /// Returns the conjugate transpose of this matrix.
         /// </summary>
-        /// <example><pre>
-        /// 1, 2, 3
-        /// 4, 5, 6  will be returned as  1, 2, 3, 4, 5, 6, 7, 8, 9
-        /// 7, 8, 9
-        /// </pre></example>
-        /// <returns>An array containing the matrix's elements.</returns>
-        public T[] ToRowWiseArray()
-        {
-            return Storage.ToRowMajorArray();
-        }
-
-        #region Implemented Interfaces
-
-#if !PORTABLE
-
-        #region ICloneable
+        /// <returns>The conjugate transpose of this matrix.</returns>
+        public abstract Matrix<T> ConjugateTranspose();
 
         /// <summary>
-        /// Creates a new object that is a copy of the current instance.
+        /// Permute the rows of a matrix according to a permutation.
         /// </summary>
-        /// <returns>
-        /// A new object that is a copy of this instance.
-        /// </returns>
-        object ICloneable.Clone()
+        /// <param name="p">The row permutation to apply to this matrix.</param>
+        public virtual void PermuteRows(Permutation p)
         {
-            return Clone();
-        }
-
-        #endregion
+            if (p.Dimension != RowCount)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "p");
+            }
 
-#endif
+            // Get a sequence of inversions from the permutation.
+            var inv = p.ToInversions();
 
-        #region IEquatable<Matrix<T>>
+            for (var i = 0; i < inv.Length; i++)
+            {
+                if (inv[i] != i)
+                {
+                    var q = inv[i];
+                    for (var j = 0; j < ColumnCount; j++)
+                    {
+                        var temp = At(q, j);
+                        At(q, j, At(i, j));
+                        At(i, j, temp);
+                    }
+                }
+            }
+        }
 
         /// <summary>
-        /// Indicates whether the current object is equal to another object of the same type.
+        /// Permute the columns of a matrix according to a permutation.
         /// </summary>
-        /// <param name="other">
-        /// An object to compare with this object.
-        /// </param>
-        /// <returns>
-        /// <c>true</c> if the current object is equal to the <paramref name="other"/> parameter; otherwise, <c>false</c>.
-        /// </returns>
-        public bool Equals(Matrix<T> other)
+        /// <param name="p">The column permutation to apply to this matrix.</param>
+        public virtual void PermuteColumns(Permutation p)
         {
-            if (other == null)
+            if (p.Dimension != ColumnCount)
             {
-                return false;
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "p");
             }
 
-            return Storage.Equals(other.Storage);
-        }
-
-        #endregion
-
-        #region IFormattable
-
-        /// <summary>
-        /// Returns a <see cref="System.String"/> that represents this instance.
-        /// </summary>
-        /// <param name="format">
-        /// The format to use.
-        /// </param>
-        /// <param name="formatProvider">
-        /// The format provider to use.
-        /// </param>
-        /// <returns>
-        /// A <see cref="System.String"/> that represents this instance.
-        /// </returns>
-        public virtual string ToString(string format, IFormatProvider formatProvider = null)
-        {
-            var stringBuilder = new StringBuilder();
-            for (var row = 0; row < RowCount; row++)
-            {
-                for (var column = 0; column < ColumnCount; column++)
-                {
-                    stringBuilder.Append(At(row, column).ToString(format, formatProvider));
-                    if (column != ColumnCount - 1)
-                    {
-                        stringBuilder.Append(formatProvider.GetTextInfo().ListSeparator);
-                    }
-                }
-
-                if (row != RowCount - 1)
-                {
-                    stringBuilder.Append(Environment.NewLine);
-                }
-            }
-
-            return stringBuilder.ToString();
-        }
-
-        #endregion
-
-        #endregion
-
-        #region System.Object overrides
-
-        /// <summary>
-        /// Determines whether the specified <see cref="System.Object"/> is equal to this instance.
-        /// </summary>
-        /// <param name="obj">The <see cref="System.Object"/> to compare with this instance.</param>
-        /// <returns>
-        ///     <c>true</c> if the specified <see cref="System.Object"/> is equal to this instance; otherwise, <c>false</c>.
-        /// </returns>
-        public override bool Equals(object obj)
-        {
-            var other = obj as Matrix<T>;
-            return other != null && Storage.Equals(other.Storage);
-        }
-
-        /// <summary>
-        /// Returns a hash code for this instance.
-        /// </summary>
-        /// <returns>
-        /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table. 
-        /// </returns>
-        public override int GetHashCode()
-        {
-            return Storage.GetHashCode();
-        }
-
-        #endregion
-
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>        
-        /// <returns>The transpose of this matrix.</returns>
-        public virtual Matrix<T> Transpose()
-        {
-            var ret = CreateMatrix(ColumnCount, RowCount);
-            for (var j = 0; j < ColumnCount; j++)
-            {
-                for (var i = 0; i < RowCount; i++)
-                {
-                    ret.At(j, i, At(i, j));
-                }
-            }
-
-            return ret;
-        }
-
-        /// <summary>
-        /// Returns the conjugate transpose of this matrix.
-        /// </summary>        
-        /// <returns>The conjugate transpose of this matrix.</returns>
-        public abstract Matrix<T> ConjugateTranspose();
-
-        /// <summary>
-        /// Permute the rows of a matrix according to a permutation.
-        /// </summary>
-        /// <param name="p">The row permutation to apply to this matrix.</param>
-        public virtual void PermuteRows(Permutation p)
-        {
-            if (p.Dimension != RowCount)
-            {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "p");
-            }
-
-            // Get a sequence of inversions from the permutation.
-            var inv = p.ToInversions();
-
-            for (var i = 0; i < inv.Length; i++)
-            {
-                if (inv[i] != i)
-                {
-                    var q = inv[i];
-                    for (var j = 0; j < ColumnCount; j++)
-                    {
-                        var temp = At(q, j);
-                        At(q, j, At(i, j));
-                        At(i, j, temp);
-                    }
-                }
-            }
-        }
-
-        /// <summary>
-        /// Permute the columns of a matrix according to a permutation.
-        /// </summary>
-        /// <param name="p">The column permutation to apply to this matrix.</param>
-        public virtual void PermuteColumns(Permutation p)
-        {
-            if (p.Dimension != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "p");
-            }
-
-            // Get a sequence of inversions from the permutation.
-            var inv = p.ToInversions();
+            // Get a sequence of inversions from the permutation.
+            var inv = p.ToInversions();
 
             for (var i = 0; i < inv.Length; i++)
             {
@@ -1501,7 +1252,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         }
 
         /// <summary>
-        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix, 
+        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix,
         /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
         /// The values of off the off diagonal matrices/blocks are set to zero.
         /// </summary>
@@ -1555,8 +1306,8 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         public abstract T L1Norm();
 
         /// <summary>Calculates the L2 norm.</summary>
-        /// <returns>The L2 norm of the matrix.</returns>   
-        /// <remarks>For sparse matrices, the L2 norm is computed using a dense implementation of singular value decomposition. 
+        /// <returns>The L2 norm of the matrix.</returns>
+        /// <remarks>For sparse matrices, the L2 norm is computed using a dense implementation of singular value decomposition.
         /// In a later release, it will be replaced with a sparse implementation.</remarks>
         public virtual T L2Norm()
         {
@@ -1568,24 +1319,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         public abstract T FrobeniusNorm();
 
         /// <summary>Calculates the infinity norm of this matrix.</summary>
-        /// <returns>The infinity norm of this matrix.</returns>   
+        /// <returns>The infinity norm of this matrix.</returns>
         public abstract T InfinityNorm();
 
-        /// <summary>
-        /// Iterates throw each element in the matrix (row-wise).
-        /// </summary>
-        /// <returns>The value at the current iteration along with its position (row, column, value).</returns>
-        public virtual IEnumerable<Tuple<int, int, T>> IndexedEnumerator()
-        {
-            for (var row = 0; row < RowCount; row++)
-            {
-                for (var column = 0; column < ColumnCount; column++)
-                {
-                    yield return new Tuple<int, int, T>(row, column, At(row, column));
-                }
-            }
-        }
-
         /// <summary>
         /// Gets a value indicating whether this matrix is symmetric.
         /// </summary>
@@ -1612,5 +1348,240 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
                 return true;
             }
         }
+
+        /// <summary>
+        /// Returns an <see cref="IEnumerator{T}"/> that enumerates over the matrix columns.
+        /// </summary>
+        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over the matrix columns</returns>
+        /// <seealso cref="IEnumerator{T}"/>
+        public virtual IEnumerable<Tuple<int, Vector<T>>> ColumnEnumerator()
+        {
+            for (var i = 0; i < ColumnCount; i++)
+            {
+                yield return new Tuple<int, Vector<T>>(i, Column(i));
+            }
+        }
+
+        /// <summary>
+        /// Returns an <see cref="IEnumerator{T}"/> that enumerates the requested matrix columns.
+        /// </summary>
+        /// <param name="index">The column to start enumerating over.</param>
+        /// <param name="length">The number of columns to enumerating over.</param>
+        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over requested matrix columns.</returns>
+        /// <seealso cref="IEnumerator{T}"/>
+        /// <exception cref="ArgumentOutOfRangeException">If:
+        /// <list><item><paramref name="index"/> is negative,
+        /// or greater than or equal to the number of columns.</item>
+        /// <item><c>(index + length) &gt;= Columns.</c></item></list>
+        /// </exception>
+        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
+        public virtual IEnumerable<Tuple<int, Vector<T>>> ColumnEnumerator(int index, int length)
+        {
+            if (index >= ColumnCount || index < 0)
+            {
+                throw new ArgumentOutOfRangeException("index");
+            }
+
+            if (index + length > ColumnCount)
+            {
+                throw new ArgumentOutOfRangeException("length");
+            }
+
+            if (length < 1)
+            {
+                throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
+            }
+
+            var maxIndex = index + length;
+            for (var i = index; i < maxIndex; i++)
+            {
+                yield return new Tuple<int, Vector<T>>(i, Column(i));
+            }
+        }
+
+        /// <summary>
+        /// Returns an <see cref="IEnumerator{T}"/> that enumerates the requested matrix rows.
+        /// </summary>
+        /// <param name="index">The row to start enumerating over.</param>
+        /// <param name="length">The number of rows to enumerating over.</param>
+        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over requested matrix rows.</returns>
+        /// <seealso cref="IEnumerator{T}"/>
+        /// <exception cref="ArgumentOutOfRangeException">If:
+        /// <list><item><paramref name="index"/> is negative,
+        /// or greater than or equal to the number of rows.</item>
+        /// <item><c>(index + length) &gt;= Rows.</c></item></list></exception>
+        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
+        public virtual IEnumerable<Tuple<int, Vector<T>>> RowEnumerator(int index, int length)
+        {
+            if (index >= RowCount || index < 0)
+            {
+                throw new ArgumentOutOfRangeException("index");
+            }
+
+            if (index + length > RowCount)
+            {
+                throw new ArgumentOutOfRangeException("length");
+            }
+
+            if (length < 1)
+            {
+                throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
+            }
+
+            var maxi = index + length;
+            for (var i = index; i < maxi; i++)
+            {
+                yield return new Tuple<int, Vector<T>>(i, Row(i));
+            }
+        }
+
+        /// <summary>
+        /// Returns an <see cref="IEnumerator{T}"/> that enumerates over the matrix rows.
+        /// </summary>
+        /// <returns>An <see cref="IEnumerator{T}"/> that enumerates over the matrix rows</returns>
+        /// <seealso cref="IEnumerator{T}"/>
+        public virtual IEnumerable<Tuple<int, Vector<T>>> RowEnumerator()
+        {
+            for (var i = 0; i < RowCount; i++)
+            {
+                yield return new Tuple<int, Vector<T>>(i, Row(i));
+            }
+        }
+
+        /// <summary>
+        /// Iterates throw each element in the matrix (row-wise).
+        /// </summary>
+        /// <returns>The value at the current iteration along with its position (row, column, value).</returns>
+        public virtual IEnumerable<Tuple<int, int, T>> IndexedEnumerator()
+        {
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    yield return new Tuple<int, int, T>(row, column, At(row, column));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Returns this matrix as a multidimensional array.
+        /// </summary>
+        /// <returns>A multidimensional containing the values of this matrix.</returns>
+        public T[,] ToArray()
+        {
+            return Storage.ToArray();
+        }
+
+        /// <summary>
+        /// Returns the matrix's elements as an array with the data laid out column-wise.
+        /// </summary>
+        /// <example><pre>
+        /// 1, 2, 3
+        /// 4, 5, 6  will be returned as  1, 4, 7, 2, 5, 8, 3, 6, 9
+        /// 7, 8, 9
+        /// </pre></example>
+        /// <returns>An array containing the matrix's elements.</returns>
+        public T[] ToColumnWiseArray()
+        {
+            return Storage.ToColumnMajorArray();
+        }
+
+        /// <summary>
+        /// Returns the matrix's elements as an array with the data laid row-wise.
+        /// </summary>
+        /// <example><pre>
+        /// 1, 2, 3
+        /// 4, 5, 6  will be returned as  1, 2, 3, 4, 5, 6, 7, 8, 9
+        /// 7, 8, 9
+        /// </pre></example>
+        /// <returns>An array containing the matrix's elements.</returns>
+        public T[] ToRowWiseArray()
+        {
+            return Storage.ToRowMajorArray();
+        }
+
+        /// <summary>
+        /// Returns a <see cref="System.String"/> that represents this instance.
+        /// </summary>
+        /// <returns>
+        /// A <see cref="System.String"/> that represents this instance.
+        /// </returns>
+        public override string ToString()
+        {
+            return ToString(null);
+        }
+
+        /// <summary>
+        /// Returns a <see cref="System.String"/> that represents this instance.
+        /// </summary>
+        /// <param name="format">
+        /// The format to use.
+        /// </param>
+        /// <param name="formatProvider">
+        /// The format provider to use.
+        /// </param>
+        /// <returns>
+        /// A <see cref="System.String"/> that represents this instance.
+        /// </returns>
+        public virtual string ToString(string format, IFormatProvider formatProvider = null)
+        {
+            var stringBuilder = new StringBuilder();
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    stringBuilder.Append(At(row, column).ToString(format, formatProvider));
+                    if (column != ColumnCount - 1)
+                    {
+                        stringBuilder.Append(formatProvider.GetTextInfo().ListSeparator);
+                    }
+                }
+
+                if (row != RowCount - 1)
+                {
+                    stringBuilder.Append(Environment.NewLine);
+                }
+            }
+
+            return stringBuilder.ToString();
+        }
+
+        /// <summary>
+        /// Returns a hash code for this instance.
+        /// </summary>
+        /// <returns>
+        /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table.
+        /// </returns>
+        public override int GetHashCode()
+        {
+            return Storage.GetHashCode();
+        }
+
+        /// <summary>
+        /// Indicates whether the current object is equal to another object of the same type.
+        /// </summary>
+        /// <param name="other">
+        /// An object to compare with this object.
+        /// </param>
+        /// <returns>
+        /// <c>true</c> if the current object is equal to the <paramref name="other"/> parameter; otherwise, <c>false</c>.
+        /// </returns>
+        public bool Equals(Matrix<T> other)
+        {
+            return other != null && Storage.Equals(other.Storage);
+        }
+
+        /// <summary>
+        /// Determines whether the specified <see cref="System.Object"/> is equal to this instance.
+        /// </summary>
+        /// <param name="obj">The <see cref="System.Object"/> to compare with this instance.</param>
+        /// <returns>
+        ///     <c>true</c> if the specified <see cref="System.Object"/> is equal to this instance; otherwise, <c>false</c>.
+        /// </returns>
+        public override bool Equals(object obj)
+        {
+            var other = obj as Matrix<T>;
+            return other != null && Storage.Equals(other.Storage);
+        }
     }
 }