LA: Clean up and extend matrix enumeration

src/Examples/LinearAlgebra/MatrixDataAccessor.cs

@@ -137,7 +137,7 @@ namespace Examples.LinearAlgebraExamples
 
             // 9. Get columns using column enumerator. If you need all columns you may use ColumnEnumerator without parameters
             Console.WriteLine(@"9. Get columns using column enumerator");
-            foreach (var keyValuePair in matrix.ColumnEnumerator(2, 4))
+            foreach (var keyValuePair in matrix.EnumerateColumnsIndexed(2, 4))
             {
                 Console.WriteLine(@"Column {0}: {1}", keyValuePair.Item1, keyValuePair.Item2.ToString("#0.00\t", formatProvider));
             }
@@ -146,7 +146,7 @@ namespace Examples.LinearAlgebraExamples
 
             // 10. Get rows using row enumerator. If you need all rows you may use RowEnumerator without parameters
             Console.WriteLine(@"10. Get rows using row enumerator");
-            foreach (var keyValuePair in matrix.RowEnumerator(4, 3))
+            foreach (var keyValuePair in matrix.EnumerateRowsIndexed(4, 3))
             {
                 Console.WriteLine(@"Row {0}: {1}", keyValuePair.Item1, keyValuePair.Item2.ToString("#0.00\t", formatProvider));
             }

src/Examples/LinearAlgebra/MatrixNorms.cs

@@ -95,7 +95,7 @@ namespace Examples.LinearAlgebraExamples
 
             // 5. Normalize matrix columns
             Console.WriteLine(@"5. Normalize matrix columns: before normalize");
-            foreach (var keyValuePair in matrix.ColumnEnumerator())
+            foreach (var keyValuePair in matrix.EnumerateColumnsIndexed())
             {
                 Console.WriteLine(@"Column {0} 2-nd norm is: {1}", keyValuePair.Item1, keyValuePair.Item2.L2Norm());
             }
@@ -103,7 +103,7 @@ namespace Examples.LinearAlgebraExamples
             Console.WriteLine();
             var normalized = matrix.NormalizeColumns(2);
             Console.WriteLine(@"5. Normalize matrix columns: after normalize");
-            foreach (var keyValuePair in normalized.ColumnEnumerator())
+            foreach (var keyValuePair in normalized.EnumerateColumnsIndexed())
             {
                 Console.WriteLine(@"Column {0} 2-nd norm is: {1}", keyValuePair.Item1, keyValuePair.Item2.L2Norm());
             }
@@ -112,7 +112,7 @@ namespace Examples.LinearAlgebraExamples
 
             // 6. Normalize matrix columns
             Console.WriteLine(@"6. Normalize matrix rows: before normalize");
-            foreach (var keyValuePair in matrix.RowEnumerator())
+            foreach (var keyValuePair in matrix.EnumerateRowsIndexed())
             {
                 Console.WriteLine(@"Row {0} 2-nd norm is: {1}", keyValuePair.Item1, keyValuePair.Item2.L2Norm());
             }
@@ -120,7 +120,7 @@ namespace Examples.LinearAlgebraExamples
             Console.WriteLine();
             normalized = matrix.NormalizeRows(2);
             Console.WriteLine(@"6. Normalize matrix rows: after normalize");
-            foreach (var keyValuePair in normalized.RowEnumerator())
+            foreach (var keyValuePair in normalized.EnumerateRowsIndexed())
             {
                 Console.WriteLine(@"Row {0} 2-nd norm is: {1}", keyValuePair.Item1, keyValuePair.Item2.L2Norm());
             }

src/FSharp/LinearAlgebra.Double.Matrix.fs

@@ -36,12 +36,6 @@ open MathNet.Numerics.LinearAlgebra
 [<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
 module Matrix =
 
-    /// Creates a sequence that iterates the non-zero entries in the matrix.
-    let inline nonZeroEntries (A: #Matrix<_>) =
-        seq { for i in 0 .. A.RowCount-1 do
-                for j in 0 .. A.ColumnCount-1 do
-                  if A.At(i,j) <> 0.0 then yield (i, j, A.At(i,j)) }
-
     /// Returns the sum of all elements of a matrix.
     let inline sum (A: #Matrix<_>) =
         let mutable f = 0.0
@@ -50,14 +44,6 @@ module Matrix =
                 f <- f + A.At(i,j)
         f
 
-    /// Returns the sum of the results generated by applying a position dependent function to each column of the matrix.
-    let inline sumColsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<_>) =
-        A.ColumnEnumerator() |> Seq.map (fun (j,col) -> f j col) |> Seq.reduce (+)
-
-    /// Returns the sum of the results generated by applying a position dependent function to each row of the matrix.
-    let inline sumRowsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<_>) =
-        A.RowEnumerator() |> Seq.map (fun (i,row) -> f i row)  |> Seq.reduce (+)
-
     /// Fold all columns into one row vector.
     let inline foldByCol (f: float -> float -> float) acc (A: #Matrix<float>) =
         let v = new DenseVector(A.ColumnCount)

src/FSharp/LinearAlgebra.Matrix.fs

@@ -266,3 +266,14 @@ module Matrix =
         for i=0 to A.ColumnCount-1 do
             macc <- f macc (A.Item(k,i))
         macc
+
+    /// Returns the sum of the results generated by applying a position dependent function to each column of the matrix.
+    let inline sumColsBy f (A: #Matrix<_>) =
+        A.EnumerateColumnsIndexed() |> Seq.map (fun (j,col) -> f j col) |> Seq.reduce (+)
+
+    /// Returns the sum of the results generated by applying a position dependent function to each row of the matrix.
+    let inline sumRowsBy f (A: #Matrix<_>) =
+        A.EnumerateRowsIndexed() |> Seq.map (fun (i,row) -> f i row)  |> Seq.reduce (+)
+
+    /// Creates a sequence that iterates the non-zero entries in the matrix.
+    let nonZeroEntries (A: #Matrix<_>) = A.EnumerateNonZero()

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -1070,44 +1070,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
-        /// <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 override IEnumerable<Tuple<int, int, Complex>> IndexedEnumerator()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (var row = 0; row < RowCount - 1; row++)
-            {
-                var start = rowPointers[row];
-                var end = rowPointers[row + 1];
-
-                if (start == end)
-                {
-                    continue;
-                }
-
-                for (var index = start; index < end; index++)
-                {
-                    yield return new Tuple<int, int, Complex>(row, columnIndices[index], values[index]);
-                }
-            }
-
-            var lastRow = rowPointers.Length - 1;
-
-            if (rowPointers[lastRow] < valueCount)
-            {
-                for (var index = rowPointers[lastRow]; index < valueCount; index++)
-                {
-                    yield return new Tuple<int, int, Complex>(lastRow, columnIndices[index], values[index]);
-                }
-            }
-        }
-
         /// <summary>
         /// Gets a value indicating whether this matrix is symmetric.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -1064,44 +1064,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
-        /// <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 override IEnumerable<Tuple<int, int, Complex32>> IndexedEnumerator()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (var row = 0; row < RowCount - 1; row++)
-            {
-                var start = rowPointers[row];
-                var end = rowPointers[row + 1];
-
-                if (start == end)
-                {
-                    continue;
-                }
-
-                for (var index = start; index < end; index++)
-                {
-                    yield return new Tuple<int, int, Complex32>(row, columnIndices[index], values[index]);
-                }
-            }
-
-            var lastRow = rowPointers.Length - 1;
-
-            if (rowPointers[lastRow] < valueCount)
-            {
-                for (var index = rowPointers[lastRow]; index < valueCount; index++)
-                {
-                    yield return new Tuple<int, int, Complex32>(lastRow, columnIndices[index], values[index]);
-                }
-            }
-        }
-
         /// <summary>
         /// Gets a value indicating whether this matrix is symmetric.
         /// </summary>

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -1088,44 +1088,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
-        /// <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 override IEnumerable<Tuple<int, int, double>> IndexedEnumerator()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (var row = 0; row < RowCount - 1; row++)
-            {
-                var start = rowPointers[row];
-                var end = rowPointers[row + 1];
-
-                if (start == end)
-                {
-                    continue;
-                }
-
-                for (var index = start; index < end; index++)
-                {
-                    yield return new Tuple<int, int, double>(row, columnIndices[index], values[index]);
-                }
-            }
-
-            var lastRow = rowPointers.Length - 1;
-
-            if (rowPointers[lastRow] < valueCount)
-            {
-                for (var index = rowPointers[lastRow]; index < valueCount; index++)
-                {
-                    yield return new Tuple<int, int, double>(lastRow, columnIndices[index], values[index]);
-                }
-            }
-        }
-
         /// <summary>
         /// Gets a value indicating whether this matrix is symmetric.
         /// </summary>

src/Numerics/LinearAlgebra/Matrix.cs

@@ -1254,97 +1254,134 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates over the matrix columns.
+        /// Returns an IEnumerable that can be used to iterate through all values of the matrix.
         /// </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()
+        /// <remarks>
+        /// The enumerator will include all values, even if they are zero.
+        /// The ordering of the values is unspecified (not necessarily column by column or row by row).
+        /// </remarks>
+        public IEnumerable<T> Enumerate()
+        {
+            return Storage.Enumerate();
+        }
+
+        /// <summary>
+        /// Returns an IEnumerable that can be used to iterate through all values of the matrix and their index.
+        /// </summary>
+        /// <remarks>
+        /// The enumerator returns a Tuple with the first two values being the row and column index
+        /// and the third value being the value of the element at that index.
+        /// The enumerator will include all values, even if they are zero.
+        /// </remarks>
+        public IEnumerable<Tuple<int, int, T>> EnumerateIndexed()
+        {
+            return Storage.EnumerateIndexed();
+        }
+
+        /// <summary>
+        /// Returns an IEnumerable that can be used to iterate through all non-zero values of the matrix and their index.
+        /// </summary>
+        /// <remarks>
+        /// The enumerator returns a Tuple with the first two values being the row and column index
+        /// and the third value being the value of the element at that index.
+        /// The enumerator will skip all elements with a zero value.
+        /// </remarks>
+        public IEnumerable<Tuple<int, int, T>> EnumerateNonZero()
+        {
+            return Storage.EnumerateNonZero();
+        }
+
+        /// <summary>
+        /// Returns an IEnumerable that can be used to iterate through all columns of the matrix.
+        /// </summary>
+        public IEnumerable<Vector<T>> EnumerateColumns()
         {
             for (var i = 0; i < ColumnCount; i++)
             {
-                yield return new Tuple<int, Vector<T>>(i, Column(i));
+                yield return Column(i);
             }
         }
 
         /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates the requested matrix columns.
+        /// Returns an IEnumerable that can be used to iterate through a subset of all columns of the matrix.
         /// </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)
+        public IEnumerable<Vector<T>> EnumerateColumns(int index, int length)
         {
-            if (index >= ColumnCount || index < 0)
+            var maxIndex = Math.Min(index + length, ColumnCount);
+            for (var i = Math.Max(index, 0); i < maxIndex; i++)
             {
-                throw new ArgumentOutOfRangeException("index");
+                yield return Column(i);
             }
+        }
 
-            if (index + length > ColumnCount)
+        /// <summary>
+        /// Returns an IEnumerable that can be used to iterate through all columns of the matrix and their index.
+        /// </summary>
+        /// <remarks>
+        /// The enumerator returns a Tuple with the first value being the column index
+        /// and the second value being the value of the column at that index.
+        /// </remarks>
+        public IEnumerable<Tuple<int, Vector<T>>> EnumerateColumnsIndexed()
+        {
+            for (var i = 0; i < ColumnCount; i++)
             {
-                throw new ArgumentOutOfRangeException("length");
+                yield return new Tuple<int, Vector<T>>(i, Column(i));
             }
+        }
 
-            if (length < 1)
+        /// <summary>
+        /// Returns an IEnumerable that can be used to iterate through a subset of all columns of the matrix and their index.
+        /// </summary>
+        /// <param name="index">The column to start enumerating over.</param>
+        /// <param name="length">The number of columns to enumerating over.</param>
+        /// <remarks>
+        /// The enumerator returns a Tuple with the first value being the column index
+        /// and the second value being the value of the column at that index.
+        /// </remarks>
+        public IEnumerable<Tuple<int, Vector<T>>> EnumerateColumnsIndexed(int index, int length)
+        {
+            var maxIndex = Math.Min(index + length, ColumnCount);
+            for (var i = Math.Max(index, 0); i < maxIndex; i++)
             {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "length");
+                yield return new Tuple<int, Vector<T>>(i, Column(i));
             }
+        }
 
-            var maxIndex = index + length;
-            for (var i = index; i < maxIndex; i++)
+        /// <summary>
+        /// Returns an IEnumerable that can be used to iterate through all rows of the matrix.
+        /// </summary>
+        public IEnumerable<Vector<T>> EnumerateRows()
+        {
+            for (var i = 0; i < RowCount; i++)
             {
-                yield return new Tuple<int, Vector<T>>(i, Column(i));
+                yield return Row(i);
             }
         }
 
         /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates the requested matrix rows.
+        /// Returns an IEnumerable that can be used to iterate through a subset of all rows of the matrix.
         /// </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)
+        public IEnumerable<Vector<T>> EnumerateRows(int index, int length)
         {
-            if (index >= RowCount || index < 0)
-            {
-                throw new ArgumentOutOfRangeException("index");
-            }
-
-            if (index + length > RowCount)
+            var maxIndex = Math.Min(index + length, RowCount);
+            for (var i = Math.Max(index, 0); i < maxIndex; i++)
             {
-                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));
+                yield return Row(i);
             }
         }
 
         /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that enumerates over the matrix rows.
+        /// Returns an IEnumerable that can be used to iterate through all rows of the matrix and their index.
         /// </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()
+        /// <remarks>
+        /// The enumerator returns a Tuple with the first value being the row index
+        /// and the second value being the value of the row at that index.
+        /// </remarks>
+        public IEnumerable<Tuple<int, Vector<T>>> EnumerateRowsIndexed()
         {
             for (var i = 0; i < RowCount; i++)
             {
@@ -1353,17 +1390,20 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Iterates through each element in the matrix (row-wise).
+        /// Returns an IEnumerable that can be used to iterate through a subset of all rows of the matrix and their index.
         /// </summary>
-        /// <returns>The value at the current iteration along with its position (row, column, value).</returns>
-        public virtual IEnumerable<Tuple<int, int, T>> IndexedEnumerator()
+        /// <param name="index">The row to start enumerating over.</param>
+        /// <param name="length">The number of rows to enumerating over.</param>
+        /// <remarks>
+        /// The enumerator returns a Tuple with the first value being the row index
+        /// and the second value being the value of the row at that index.
+        /// </remarks>
+        public IEnumerable<Tuple<int, Vector<T>>> EnumerateRowsIndexed(int index, int length)
         {
-            for (var row = 0; row < RowCount; row++)
+            var maxIndex = Math.Min(index + length, RowCount);
+            for (var i = Math.Max(index, 0); i < maxIndex; i++)
             {
-                for (var column = 0; column < ColumnCount; column++)
-                {
-                    yield return new Tuple<int, int, T>(row, column, At(row, column));
-                }
+                yield return new Tuple<int, Vector<T>>(i, Row(i));
             }
         }
 

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -1087,44 +1087,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
-        /// <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 override IEnumerable<Tuple<int, int, float>> IndexedEnumerator()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-            var valueCount = _storage.ValueCount;
-
-            for (var row = 0; row < RowCount - 1; row++)
-            {
-                var start = rowPointers[row];
-                var end = rowPointers[row + 1];
-
-                if (start == end)
-                {
-                    continue;
-                }
-
-                for (var index = start; index < end; index++)
-                {
-                    yield return new Tuple<int, int, float>(row, columnIndices[index], values[index]);
-                }
-            }
-
-            var lastRow = rowPointers.Length - 1;
-
-            if (rowPointers[lastRow] < valueCount)
-            {
-                for (var index = rowPointers[lastRow]; index < valueCount; index++)
-                {
-                    yield return new Tuple<int, int, float>(lastRow, columnIndices[index], values[index]);
-                }
-            }
-        }
-
         /// <summary>
         /// Gets a value indicating whether this matrix is symmetric.
         /// </summary>

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -30,6 +30,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Threading;
 
@@ -276,6 +277,43 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return new DenseColumnMajorMatrixStorage<T>(rows, columns, array);
         }
 
+        // ENUMERATION
+
+        public override IEnumerable<T> Enumerate()
+        {
+            return Data;
+        }
+
+        public override IEnumerable<Tuple<int, int, T>> EnumerateIndexed()
+        {
+            int index = 0;
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    yield return new Tuple<int, int, T>(i, j, Data[index]);
+                    index++;
+                }
+            }
+        }
+
+        public override IEnumerable<Tuple<int, int, T>> EnumerateNonZero()
+        {
+            int index = 0;
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    var x = Data[index];
+                    if (!Zero.Equals(x))
+                    {
+                        yield return new Tuple<int, int, T>(i, j, x);
+                    }
+                    index++;
+                }
+            }
+        }
+
         // MATRIX COPY
 
         internal override void CopyToUnchecked(MatrixStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@@ -244,6 +244,45 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        // ENUMERATION
+
+        public override IEnumerable<T> Enumerate()
+        {
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    // PERF: consider to break up loop to avoid branching
+                    yield return i == j ? Data[i] : Zero;
+                }
+            }
+        }
+
+        public override IEnumerable<Tuple<int, int, T>> EnumerateIndexed()
+        {
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    // PERF: consider to break up loop to avoid branching
+                    yield return i == j
+                        ? new Tuple<int, int, T>(i, i, Data[i])
+                        : new Tuple<int, int, T>(i, j, Zero);
+                }
+            }
+        }
+
+        public override IEnumerable<Tuple<int, int, T>> EnumerateNonZero()
+        {
+            for (int i = 0; i < Data.Length; i++)
+            {
+                if (!Zero.Equals(Data[i]))
+                {
+                    yield return new Tuple<int, int, T>(i, i, Data[i]);
+                }
+            }
+        }
+
         // MATRIX COPY
 
         internal override void CopyToUnchecked(MatrixStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 using System;
+using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Storage
@@ -226,6 +227,45 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return hash;
         }
 
+        // ENUMERATION
+
+        public virtual IEnumerable<T> Enumerate()
+        {
+            for (int i = 0; i < RowCount; i++)
+            {
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    yield return At(i, j);
+                }
+            }
+        }
+
+        public virtual IEnumerable<Tuple<int, int, T>> EnumerateIndexed()
+        {
+            for (int i = 0; i < RowCount; i++)
+            {
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    yield return new Tuple<int, int, T>(i, j, At(i, j));
+                }
+            }
+        }
+
+        public virtual IEnumerable<Tuple<int, int, T>> EnumerateNonZero()
+        {
+            for (int i = 0; i < RowCount; i++)
+            {
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    var x = At(i, j);
+                    if (!Zero.Equals(x))
+                    {
+                        yield return new Tuple<int, int, T>(i, j, x);
+                    }
+                }
+            }
+        }
+
         // MATRIX COPY
 
         public void CopyTo(MatrixStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -638,6 +638,52 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        // ENUMERATION
+
+        public override IEnumerable<T> Enumerate()
+        {
+            int k = 0;
+            for (int row = 0; row < RowCount; row++)
+            {
+                for (int col = 0; col < ColumnCount; col++)
+                {
+                    yield return k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (ColumnIndices[k]) == col
+                        ? Values[k++]
+                        : Zero;
+                }
+            }
+        }
+
+        public override IEnumerable<Tuple<int, int, T>> EnumerateIndexed()
+        {
+            int k = 0;
+            for (int row = 0; row < RowCount; row++)
+            {
+                for (int col = 0; col < ColumnCount; col++)
+                {
+                    yield return k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (ColumnIndices[k]) == col
+                        ? new Tuple<int, int, T>(row, col, Values[k++])
+                        : new Tuple<int, int, T>(row, col, Zero);
+                }
+            }
+        }
+
+        public override IEnumerable<Tuple<int, int, T>> EnumerateNonZero()
+        {
+            for (int row = 0; row < RowCount; row++)
+            {
+                var startIndex = RowPointers[row];
+                var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
+                for (var j = startIndex; j < endIndex; j++)
+                {
+                    if (!Zero.Equals(Values[j]))
+                    {
+                        yield return new Tuple<int, int, T>(row, ColumnIndices[j], Values[j]);
+                    }
+                }
+            }
+        }
+
         // MATRIX COPY
 
         internal override void CopyToUnchecked(MatrixStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Vector.cs

@@ -303,7 +303,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// Returns an IEnumerable that can be used to iterate through all values of the vector and their index.
         /// </summary>
         /// <remarks>
-        /// The enumerator returns a <seealso cref="Tuple{T,K}"/> with the first value being the element index
+        /// The enumerator returns a Tuple with the first value being the element index
         /// and the second value being the value of the element at that index.
         /// The enumerator will include all values, even if they are zero.
         /// </remarks>
@@ -316,7 +316,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// Returns an IEnumerable that can be used to iterate through all non-zero values of the vector and their index.
         /// </summary>
         /// <remarks>
-        /// The enumerator returns a <seealso cref="Tuple{T,K}"/> with the first value being the element index
+        /// The enumerator returns a Tuple with the first value being the element index
         /// and the second value being the value of the element at that index.
         /// The enumerator will skip all elements with a zero value.
         /// </remarks>

src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.cs

@@ -365,17 +365,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests
         public void CanEnumerateWithIndex()
         {
             var dense = CreateDenseRandom(2, 3, 0);
-            using (var enumerator = dense.IndexedEnumerator().GetEnumerator())
-                for (int i = 0; i < 2; i++)
-                {
-                    for (int j = 0; j < 3; j++)
-                    {
-                        enumerator.MoveNext();
-                        Assert.AreEqual(i, enumerator.Current.Item1);
-                        Assert.AreEqual(j, enumerator.Current.Item2);
-                        Assert.AreEqual(dense[i, j], enumerator.Current.Item3);
-                    }
-                }
+            int rowIdxSum = 0, colIdxSum = 0;
+            foreach (var value in dense.EnumerateIndexed())
+            {
+                rowIdxSum += value.Item1;
+                colIdxSum += value.Item2;
+                Assert.AreEqual(dense[value.Item1, value.Item2], value.Item3);
+            }
+            Assert.AreEqual(dense.RowCount*(dense.RowCount - 1)/2*dense.ColumnCount, rowIdxSum);
+            Assert.AreEqual(dense.ColumnCount*(dense.ColumnCount - 1)/2*dense.RowCount, colIdxSum);
         }
     }
 }