LA: Vector enumeration migrated to storage, unified Zero definition

src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// If <paramref name="order"/> is less than one.
         /// </exception>
         public DiagonalMatrix(int order)
-            : this(new DiagonalMatrixStorage<Complex>(order, order, Complex.Zero))
+            : this(new DiagonalMatrixStorage<Complex>(order, order))
          {
          }
 
@@ -85,7 +85,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// The number of columns.
         /// </param>
         public DiagonalMatrix(int rows, int columns)
-            : this(new DiagonalMatrixStorage<Complex>(rows, columns, Complex.Zero))
+            : this(new DiagonalMatrixStorage<Complex>(rows, columns))
         {
         }
 
@@ -116,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <param name="columns">The number of columns.</param>
         /// <param name="diagonalArray">The one dimensional array which contain diagonal elements.</param>
         public DiagonalMatrix(int rows, int columns, Complex[] diagonalArray)
-            : this(new DiagonalMatrixStorage<Complex>(rows, columns, Complex.Zero, diagonalArray))
+            : this(new DiagonalMatrixStorage<Complex>(rows, columns, diagonalArray))
         {
         }
 

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -75,7 +75,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// The number of columns.
         /// </param>
         public SparseMatrix(int rows, int columns)
-            : this(new SparseCompressedRowMatrixStorage<Complex>(rows, columns, Complex.Zero))
+            : this(new SparseCompressedRowMatrixStorage<Complex>(rows, columns))
         {
         }
 
@@ -494,7 +494,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseCompressedRowMatrixStorage<Complex>(ColumnCount, RowCount, Complex.Zero)
+            var ret = new SparseCompressedRowMatrixStorage<Complex>(ColumnCount, RowCount)
                 {
                     ColumnIndices = new int[valueCount],
                     Values = new Complex[valueCount]
@@ -601,7 +601,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var m = new SparseCompressedRowMatrixStorage<Complex>(order, order, Complex.Zero)
+            var m = new SparseCompressedRowMatrixStorage<Complex>(order, order)
                 {
                     ValueCount = order,
                     Values = new Complex[order],

src/Numerics/LinearAlgebra/Complex/SparseVector.cs

@@ -1044,27 +1044,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
             return base.ToString(format, formatProvider);
         }
-
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that contains the position and value of the element.
-        /// </summary>
-        /// <returns>
-        /// An <see cref="IEnumerator{T}"/> over this vector that contains the position and value of each
-        /// element.
-        /// </returns>
-        /// <remarks>
-        /// The enumerator returns a 
-        /// <seealso cref="Tuple{T,K}"/>
-        /// with the first value being the element index and the second value 
-        /// being the value of the element at that index. For sparse vectors, the enumerator will exclude all elements
-        /// with a zero value.
-        /// </remarks>
-        public override IEnumerable<Tuple<int, Complex>> GetIndexedEnumerator()
-        {
-            for (var i = 0; i < _storage.ValueCount; i++)
-            {
-                yield return new Tuple<int, Complex>(_storage.Indices[i], _storage.Values[i]);
-            }
-        }
     }
 }

src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// If <paramref name="order"/> is less than one.
         /// </exception>
         public DiagonalMatrix(int order)
-            : this(new DiagonalMatrixStorage<Complex32>(order, order, Complex32.Zero))
+            : this(new DiagonalMatrixStorage<Complex32>(order, order))
         {
         }
 
@@ -85,7 +85,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// The number of columns.
         /// </param>
         public DiagonalMatrix(int rows, int columns)
-            : this(new DiagonalMatrixStorage<Complex32>(rows, columns, Complex32.Zero))
+            : this(new DiagonalMatrixStorage<Complex32>(rows, columns))
         {
         }
 
@@ -116,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <param name="columns">The number of columns.</param>
         /// <param name="diagonalArray">The one dimensional array which contain diagonal elements.</param>
         public DiagonalMatrix(int rows, int columns, Complex32[] diagonalArray)
-            : this(new DiagonalMatrixStorage<Complex32>(rows, columns, Complex32.Zero, diagonalArray))
+            : this(new DiagonalMatrixStorage<Complex32>(rows, columns, diagonalArray))
         {
         }
 

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -75,7 +75,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// The number of columns.
         /// </param>
         public SparseMatrix(int rows, int columns)
-            : this(new SparseCompressedRowMatrixStorage<Complex32>(rows, columns, Complex32.Zero))
+            : this(new SparseCompressedRowMatrixStorage<Complex32>(rows, columns))
         {
         }
 
@@ -494,7 +494,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseCompressedRowMatrixStorage<Complex32>(ColumnCount, RowCount, Complex32.Zero)
+            var ret = new SparseCompressedRowMatrixStorage<Complex32>(ColumnCount, RowCount)
                 {
                     ColumnIndices = new int[valueCount],
                     Values = new Complex32[valueCount]
@@ -601,7 +601,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var m = new SparseCompressedRowMatrixStorage<Complex32>(order, order, Complex32.Zero)
+            var m = new SparseCompressedRowMatrixStorage<Complex32>(order, order)
                 {
                     ValueCount = order,
                     Values = new Complex32[order],

src/Numerics/LinearAlgebra/Complex32/SparseVector.cs

@@ -1044,27 +1044,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
             return base.ToString(format, formatProvider);
         }
-
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that contains the position and value of the element.
-        /// </summary>
-        /// <returns>
-        /// An <see cref="IEnumerator{T}"/> over this vector that contains the position and value of each
-        /// element.
-        /// </returns>
-        /// <remarks>
-        /// The enumerator returns a 
-        /// <seealso cref="Tuple{T,K}"/>
-        /// with the first value being the element index and the second value 
-        /// being the value of the element at that index. For sparse vectors, the enumerator will exclude all elements
-        /// with a zero value.
-        /// </remarks>
-        public override IEnumerable<Tuple<int, Complex32>> GetIndexedEnumerator()
-        {
-            for (var i = 0; i < _storage.ValueCount; i++)
-            {
-                yield return new Tuple<int, Complex32>(_storage.Indices[i], _storage.Values[i]);
-            }
-        }
     }
 }

src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@@ -70,7 +70,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// If <paramref name="order"/> is less than one.
         /// </exception>
         public DiagonalMatrix(int order)
-            : this(new DiagonalMatrixStorage<double>(order, order, 0d))
+            : this(new DiagonalMatrixStorage<double>(order, order))
          {
          }
 
@@ -84,7 +84,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// The number of columns.
         /// </param>
         public DiagonalMatrix(int rows, int columns)
-            : this(new DiagonalMatrixStorage<double>(rows, columns, 0d))
+            : this(new DiagonalMatrixStorage<double>(rows, columns))
         {
         }
 
@@ -115,7 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="columns">The number of columns.</param>
         /// <param name="diagonalArray">The one dimensional array which contain diagonal elements.</param>
         public DiagonalMatrix(int rows, int columns, double[] diagonalArray)
-            : this(new DiagonalMatrixStorage<double>(rows, columns, 0d, diagonalArray))
+            : this(new DiagonalMatrixStorage<double>(rows, columns, diagonalArray))
         {
         }
 

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -74,7 +74,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// The number of columns.
         /// </param>
         public SparseMatrix(int rows, int columns)
-            : this(new SparseCompressedRowMatrixStorage<double>(rows, columns, 0d))
+            : this(new SparseCompressedRowMatrixStorage<double>(rows, columns))
         {
         }
         
@@ -493,7 +493,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseCompressedRowMatrixStorage<double>(ColumnCount, RowCount, 0d)
+            var ret = new SparseCompressedRowMatrixStorage<double>(ColumnCount, RowCount)
                 {
                     ColumnIndices = new int[valueCount],
                     Values = new double[valueCount]
@@ -599,7 +599,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var m = new SparseCompressedRowMatrixStorage<double>(order, order, 0d)
+            var m = new SparseCompressedRowMatrixStorage<double>(order, order)
                 {
                     ValueCount = order,
                     Values = new double[order],

src/Numerics/LinearAlgebra/Double/SparseVector.cs

@@ -1084,27 +1084,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             return base.ToString(format, formatProvider);
         }
-
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that contains the position and value of the element.
-        /// </summary>
-        /// <returns>
-        /// An <see cref="IEnumerator{T}"/> over this vector that contains the position and value of each
-        /// element.
-        /// </returns>
-        /// <remarks>
-        /// The enumerator returns a 
-        /// <seealso cref="Tuple{T,K}"/>
-        /// with the first value being the element index and the second value 
-        /// being the value of the element at that index. For sparse vectors, the enumerator will exclude all elements
-        /// with a zero value.
-        /// </remarks>
-        public override IEnumerable<Tuple<int, double>> GetIndexedEnumerator()
-        {
-            for (var i = 0; i < _storage.ValueCount; i++)
-            {
-                yield return new Tuple<int, double>(_storage.Indices[i], _storage.Values[i]);
-            }
-        }
     }
 }

src/Numerics/LinearAlgebra/Generic/Common.cs

@@ -29,7 +29,7 @@
 // </copyright>
 
 
-namespace MathNet.Numerics.LinearAlgebra.Generic
+namespace MathNet.Numerics.LinearAlgebra
 {
     using System;
 
@@ -43,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// </summary>
         /// <typeparam name="T">The type to return the value of 1.0 of.</typeparam>
         /// <returns>The value of <c>1.0</c> for type T.</returns>
-        public static T SetOne<T>()
+        public static T OneOf<T>()
         {
             if (typeof(T) == typeof(System.Numerics.Complex))
             {
@@ -57,7 +57,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
 
             if (typeof(T) == typeof(double))
             {
-                return (T)(object)1.0;
+                return (T)(object)1.0d;
             }
 
             if (typeof(T) == typeof(float))
@@ -67,5 +67,35 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
 
             throw new NotSupportedException();
         }
+
+        /// <summary>
+        /// Sets the value of <c>0.0</c> for type T.
+        /// </summary>
+        /// <typeparam name="T">The type to return the value of 0.0 of.</typeparam>
+        /// <returns>The value of <c>0.0</c> for type T.</returns>
+        public static T ZeroOf<T>()
+        {
+            if (typeof(T) == typeof(System.Numerics.Complex))
+            {
+                return (T)(object)System.Numerics.Complex.Zero;
+            }
+
+            if (typeof(T) == typeof(Numerics.Complex32))
+            {
+                return (T)(object)Numerics.Complex32.Zero;
+            }
+
+            if (typeof(T) == typeof(double))
+            {
+                return (T)(object)0.0d;
+            }
+
+            if (typeof(T) == typeof(float))
+            {
+                return (T)(object)0.0f;
+            }
+
+            return default(T);
+        }
     }
 }

src/Numerics/LinearAlgebra/Generic/Factorization/LU.cs

@@ -48,7 +48,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization
         /// <summary>
         /// Value of one for T.
         /// </summary>
-        private static readonly T One = Common.SetOne<T>();
+        private static readonly T One = Common.OneOf<T>();
 
         /// <summary>
         /// Gets or sets both the L and U factors in the same matrix.

src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs

@@ -32,7 +32,6 @@
 namespace MathNet.Numerics.LinearAlgebra.Generic
 {
     using System;
-    using Distributions;
     using Factorization;
     using Properties;
 
@@ -44,12 +43,12 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <summary>
         /// The value of 1.0.
         /// </summary>
-        private static readonly T One = Common.SetOne<T>();
+        static readonly T One = Common.OneOf<T>();
 
         /// <summary>
         /// The value of 0.0.
         /// </summary>
-        private static readonly T Zero = default(T);
+        static readonly T Zero = Common.ZeroOf<T>();
 
         /// <summary>
         /// Adds another matrix to this matrix.

src/Numerics/LinearAlgebra/Generic/Matrix.cs

@@ -586,7 +586,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, row >= column ? At(row, column) : default(T));
+                    result.At(row, column, row >= column ? At(row, column) : Zero);
                 }
             }
         }
@@ -613,7 +613,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, row <= column ? At(row, column) : default(T));
+                    result.At(row, column, row <= column ? At(row, column) : Zero);
                 }
             }
         }
@@ -702,7 +702,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, row > column ? At(row, column) : default(T));
+                    result.At(row, column, row > column ? At(row, column) : Zero);
                 }
             }
         }
@@ -749,7 +749,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, row < column ? At(row, column) : default(T));
+                    result.At(row, column, row < column ? At(row, column) : Zero);
                 }
             }
         }

src/Numerics/LinearAlgebra/Generic/Vector.cs

@@ -43,22 +43,18 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
     /// <typeparam name="T">Supported data types are double, single, <see cref="Complex"/>, and <see cref="Complex32"/>.</typeparam>
     [Serializable]
     public abstract partial class Vector<T> :
-#if PORTABLE
         IFormattable, IEnumerable<T>, IEquatable<Vector<T>>
-#else
-        IFormattable, IEnumerable<T>, IEquatable<Vector<T>>
-#endif
         where T : struct, IEquatable<T>, IFormattable
     {
         /// <summary>
         /// The zero value for type T.
         /// </summary>
-        private static readonly T Zero = default(T);
+        private static readonly T Zero = Common.ZeroOf<T>();
 
         /// <summary>
         /// The value of 1.0 for type T.
         /// </summary>
-        private static readonly T One = Common.SetOne<T>();
+        private static readonly T One = Common.OneOf<T>();
 
         /// <summary>
         /// Initializes a new instance of the Vector class.
@@ -1239,34 +1235,26 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// <returns>
         /// A <see cref="T:System.Collections.Generic.IEnumerator`1"/> that can be used to iterate through the collection.
         /// </returns>
-        public virtual IEnumerator<T> GetEnumerator()
+        public IEnumerator<T> GetEnumerator()
         {
-            for (var index = 0; index < Count; index++)
-            {
-                yield return this[index];
-            }
+            return Storage.Enumerate().GetEnumerator();
         }
 
         /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that contains the position and value of the element.
+        /// Returns an <see cref="IEnumerable{T}"/> that contains the position and value of the element, for all non-zero elements.
         /// </summary>
         /// <returns>
-        /// An <see cref="IEnumerator{T}"/> over this vector that contains the position and value of each
-        /// element.
+        /// An <see cref="IEnumerable{T}"/> over this vector that contains the position and value of each element.
         /// </returns>
         /// <remarks>
-        /// The enumerator returns a
-        /// <seealso cref="Tuple{T,K}"/>
+        /// The enumerator returns a <seealso cref="Tuple{T,K}"/>
         /// with the first value being the element index and the second value
-        /// being the value of the element at that index. For sparse vectors, the enumerator will exclude all elements
-        /// with a zero value.
+        /// being the value of the element at that index.
+        /// The enumerator will exclude all elements with a zero value.
         /// </remarks>
-        public virtual IEnumerable<Tuple<int, T>> GetIndexedEnumerator()
+        public IEnumerable<Tuple<int, T>> GetIndexedEnumerator()
         {
-            for (var i = 0; i < Count; i++)
-            {
-                yield return new Tuple<int, T>(i, this[i]);
-            }
+            return Storage.EnumerateNonZero();
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@@ -70,7 +70,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// If <paramref name="order"/> is less than one.
         /// </exception>
         public DiagonalMatrix(int order)
-            : this(new DiagonalMatrixStorage<float>(order, order, 0f))
+            : this(new DiagonalMatrixStorage<float>(order, order))
         {
         }
 
@@ -84,7 +84,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// The number of columns.
         /// </param>
         public DiagonalMatrix(int rows, int columns)
-            : this(new DiagonalMatrixStorage<float>(rows, columns, 0f))
+            : this(new DiagonalMatrixStorage<float>(rows, columns))
         {
         }
 
@@ -115,7 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <param name="columns">The number of columns.</param>
         /// <param name="diagonalArray">The one dimensional array which contain diagonal elements.</param>
         public DiagonalMatrix(int rows, int columns, float[] diagonalArray)
-            : this(new DiagonalMatrixStorage<float>(rows, columns, 0f, diagonalArray))
+            : this(new DiagonalMatrixStorage<float>(rows, columns, diagonalArray))
         {
         }
 

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -74,7 +74,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// The number of columns.
         /// </param>
         public SparseMatrix(int rows, int columns)
-            : this(new SparseCompressedRowMatrixStorage<float>(rows, columns, 0f))
+            : this(new SparseCompressedRowMatrixStorage<float>(rows, columns))
         {
         }
         
@@ -493,7 +493,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             var values = _storage.Values;
             var valueCount = _storage.ValueCount;
 
-            var ret = new SparseCompressedRowMatrixStorage<float>(ColumnCount, RowCount, 0f)
+            var ret = new SparseCompressedRowMatrixStorage<float>(ColumnCount, RowCount)
                 {
                     ColumnIndices = new int[valueCount],
                     Values = new float[valueCount]
@@ -599,7 +599,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </exception>
         public static SparseMatrix Identity(int order)
         {
-            var mStorage = new SparseCompressedRowMatrixStorage<float>(order, order, 0f)
+            var mStorage = new SparseCompressedRowMatrixStorage<float>(order, order)
                 {
                     ValueCount = order,
                     Values = new float[order],

src/Numerics/LinearAlgebra/Single/SparseVector.cs

@@ -1088,27 +1088,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
             return base.ToString(format, formatProvider);
         }
-
-        /// <summary>
-        /// Returns an <see cref="IEnumerator{T}"/> that contains the position and value of the element.
-        /// </summary>
-        /// <returns>
-        /// An <see cref="IEnumerator{T}"/> over this vector that contains the position and value of each
-        /// element.
-        /// </returns>
-        /// <remarks>
-        /// The enumerator returns a 
-        /// <seealso cref="Tuple{T,K}"/>
-        /// with the first value being the element index and the second value 
-        /// being the value of the element at that index. For sparse vectors, the enumerator will exclude all elements
-        /// with a zero value.
-        /// </remarks>
-        public override IEnumerable<Tuple<int, float>> GetIndexedEnumerator()
-        {
-            for (var i = 0; i < _storage.ValueCount; i++)
-            {
-                yield return new Tuple<int, float>(_storage.Indices[i], _storage.Values[i]);
-            }
-        }
     }
 }

src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 using System;
+using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Storage
@@ -89,6 +90,24 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Clear(Data, index, count);
         }
 
+        // ENUMERATION
+
+        public override IEnumerable<T> Enumerate()
+        {
+            return Data;
+        }
+
+        public override IEnumerable<Tuple<int, T>> EnumerateNonZero()
+        {
+            for (var i = 0; i < Data.Length; i++)
+            {
+                if (!Zero.Equals(Data[i]))
+                {
+                    yield return new Tuple<int, T>(i, Data[i]);
+                }
+            }
+        }
+
         // VECTOR COPY
 
         internal override void CopyToUnchecked(VectorStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@@ -40,17 +40,15 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
     {
         // [ruegg] public fields are OK here
 
-        readonly T _zero;
         public readonly T[] Data;
 
-        internal DiagonalMatrixStorage(int rows, int columns, T zero)
+        internal DiagonalMatrixStorage(int rows, int columns)
             : base(rows, columns)
         {
-            _zero = zero;
             Data = new T[Math.Min(rows, columns)];
         }
 
-        internal DiagonalMatrixStorage(int rows, int columns, T zero, T[] data)
+        internal DiagonalMatrixStorage(int rows, int columns, T[] data)
             : base(rows, columns)
         {
             if (data == null)
@@ -63,7 +61,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, Math.Min(rows, columns)));
             }
 
-            _zero = zero;
             Data = data;
         }
 
@@ -72,7 +69,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         /// </summary>
         public override T At(int row, int column)
         {
-            return row == column ? Data[row] : _zero;
+            return row == column ? Data[row] : Zero;
         }
 
         /// <summary>
@@ -84,7 +81,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 Data[row] = value;
             }
-            else if (!_zero.Equals(value))
+            else if (!Zero.Equals(value))
             {
                 throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
             }
@@ -278,7 +275,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             if (sourceRowIndex - sourceColumnIndex != targetRowIndex - targetColumnIndex)
             {
-                if (Data.Any(x => !_zero.Equals(x)))
+                if (Data.Any(x => !Zero.Equals(x)))
                 {
                     throw new NotSupportedException();
                 }

src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs

@@ -39,6 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
     {
         // [ruegg] public fields are OK here
 
+        protected static readonly T Zero = Common.ZeroOf<T>();
         public readonly int RowCount;
         public readonly int ColumnCount;
 
@@ -133,7 +134,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 for (var j = 0; j < ColumnCount; j++)
                 {
-                    At(i, j, default(T));
+                    At(i, j, Zero);
                 }
             }
         }
@@ -144,7 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 for (var j = columnIndex; j < columnIndex + columnCount; j++)
                 {
-                    At(i, j, default(T));
+                    At(i, j, Zero);
                 }
             }
         }

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -40,8 +40,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
     {
         // [ruegg] public fields are OK here
 
-        readonly T _zero;
-
         /// <summary>
         /// The array containing the row indices of the existing rows. Element "j" of the array gives the index of the 
         /// element in the <see cref="Values"/> array that is first non-zero element in a row "j"
@@ -66,10 +64,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         /// <value>The number of non zero elements.</value>
         public int ValueCount;
 
-        internal SparseCompressedRowMatrixStorage(int rows, int columns, T zero = default(T))
+        internal SparseCompressedRowMatrixStorage(int rows, int columns)
             : base(rows, columns)
         {
-            _zero = zero;
             RowPointers = new int[rows];
             ColumnIndices = new int[0];
             Values = new T[0];
@@ -92,7 +89,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         public override T At(int row, int column)
         {
             var index = FindItem(row, column);
-            return index >= 0 ? Values[index] : _zero;
+            return index >= 0 ? Values[index] : Zero;
         }
 
         /// <summary>
@@ -108,7 +105,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             if (index >= 0)
             {
                 // Non-zero item found in matrix
-                if (_zero.Equals(value))
+                if (Zero.Equals(value))
                 {
                     // Delete existing item
                     RemoveAtIndexUnchecked(index, row);
@@ -122,7 +119,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             else
             {
                 // Item not found. Add new value
-                if (_zero.Equals(value))
+                if (Zero.Equals(value))
                 {
                     return;
                 }
@@ -575,7 +572,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             for (int i = sourceColumnIndex, j = 0; i < sourceColumnIndex + columnCount; i++, j++)
             {
                 var index = FindItem(rowIndex, i);
-                target.At(j, index >= 0 ? Values[index] : _zero);
+                target.At(j, index >= 0 ? Values[index] : Zero);
             }
         }
 

src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 using System;
+using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Storage
@@ -39,8 +40,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
     {
         // [ruegg] public fields are OK here
 
-        readonly T _zero;
-
         /// <summary>
         /// Array that contains the indices of the non-zero values.
         /// </summary>
@@ -56,10 +55,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         /// </summary>
         public int ValueCount;
 
-        internal SparseVectorStorage(int length, T zero = default(T))
+        internal SparseVectorStorage(int length)
             : base(length)
         {
-            _zero = zero;
             Indices = new int[0];
             Values = new T[0];
             ValueCount = 0;
@@ -72,7 +70,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             // Search if item idex exists in NonZeroIndices array in range "0 - nonzero values count"
             var itemIndex = Array.BinarySearch(Indices, 0, ValueCount, index);
-            return itemIndex >= 0 ? Values[itemIndex] : _zero;
+            return itemIndex >= 0 ? Values[itemIndex] : Zero;
         }
 
         /// <summary>
@@ -85,7 +83,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             if (itemIndex >= 0)
             {
                 // Non-zero item found in matrix
-                if (_zero.Equals(value))
+                if (Zero.Equals(value))
                 {
                     // Delete existing item
                     RemoveAtIndexUnchecked(itemIndex);
@@ -99,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             else
             {
                 // Item not found. Add new value
-                if (!_zero.Equals(value))
+                if (!Zero.Equals(value))
                 {
                     InsertAtIndexUnchecked(~itemIndex, index, value);
                 }
@@ -235,7 +233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 if (j >= otherSparse.ValueCount || i < ValueCount && Indices[i] < otherSparse.Indices[j])
                 {
-                    if (!_zero.Equals(Values[i++]))
+                    if (!Zero.Equals(Values[i++]))
                     {
                         return false;
                     }
@@ -244,7 +242,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
                 if (i >= ValueCount || j < otherSparse.ValueCount && otherSparse.Indices[j] < Indices[i])
                 {
-                    if (!_zero.Equals(otherSparse.Values[j++]))
+                    if (!Zero.Equals(otherSparse.Values[j++]))
                     {
                         return false;
                     }
@@ -284,6 +282,30 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return hash;
         }
 
+        // ENUMERATION
+
+        public override IEnumerable<T> Enumerate()
+        {
+            int k = 0;
+            for (int i = 0; i < Length; i++)
+            {
+                yield return k < ValueCount && Indices[k] == i
+                    ? Values[k++]
+                    : Zero;
+            }
+        }
+
+        public override IEnumerable<Tuple<int, T>> EnumerateNonZero()
+        {
+            for (var i = 0; i < ValueCount; i++)
+            {
+                if (!Zero.Equals(Values[i]))
+                {
+                    yield return new Tuple<int, T>(Indices[i], Values[i]);
+                }
+            }
+        }
+
         // VECTOR COPY
 
         internal override void CopyToUnchecked(VectorStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Storage/VectorStorage.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 using System;
+using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.LinearAlgebra.Storage
@@ -39,6 +40,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
     {
         // [ruegg] public fields are OK here
 
+        protected static readonly T Zero = Common.ZeroOf<T>();
         public readonly int Length;
 
         protected VectorStorage(int length)
@@ -113,7 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             for (var i = 0; i < Length; i++)
             {
-                At(i, default(T));
+                At(i, Zero);
             }
         }
 
@@ -121,7 +123,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             for (var i = index; i < index + count; i++)
             {
-                At(i, default(T));
+                At(i, Zero);
             }
         }
 
@@ -196,6 +198,28 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return hash;
         }
 
+        // ENUMERATION
+
+        public virtual IEnumerable<T> Enumerate()
+        {
+            for (var i = 0; i < Length; i++)
+            {
+                yield return At(i);
+            }
+        }
+
+        public virtual IEnumerable<Tuple<int, T>> EnumerateNonZero()
+        {
+            for (var i = 0; i < Length; i++)
+            {
+                var x = At(i);
+                if (!Zero.Equals(x))
+                {
+                    yield return new Tuple<int, T>(i, x);
+                }
+            }
+        }
+
         // VECTOR COPY
 
         public void CopyTo(VectorStorage<T> target, bool skipClearing = false)