LA: Functional matrix init

src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@@ -138,27 +138,40 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// Create a new dense matrix as a copy of the given enumerable.
         /// The enumerable is assumed to be in column-major order (column by column).
         /// This new matrix will be independent from the enumerable.
-        /// A new memory block will be allocated for storing the vector.
+        /// A new memory block will be allocated for storing the matrix.
         /// </summary>
         public static DenseMatrix OfColumnMajor(int rows, int columns, IEnumerable<Complex> columnMajor)
         {
             return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static DenseMatrix Create(int rows, int columns, Func<int, int, Complex> init)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix with values sampled from the provided random distribution.
+        /// </summary>
+        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfInit(rows, columns,
+                (i, j) => new Complex(distribution.Sample(), distribution.Sample())));
+        }
+
         /// <summary>
         /// Create a new dense matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Scheduled for removal in v3.0.")]
+        [Obsolete("Use DenseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public DenseMatrix(int rows, int columns, Complex value)
-            : this(rows, columns)
+            : this(DenseColumnMajorMatrixStorage<Complex>.OfInit(rows, columns, (i, j) => value))
         {
-            for (var i = 0; i < _values.Length; i++)
-            {
-                _values[i] = value;
-            }
         }
 
         /// <summary>
@@ -183,19 +196,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
         }
 
-        /// <summary>
-        /// Create a new dense matrix with values sampled from the provided random distribution.
-        /// </summary>
-        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<Complex>(rows, columns);
-            for (var i = 0; i < storage.Data.Length; i++)
-            {
-                storage.Data[i] = new Complex(distribution.Sample(), distribution.Sample());
-            }
-            return new DenseMatrix(storage);
-        }
-
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>

src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@@ -30,6 +30,7 @@
 
 namespace MathNet.Numerics.LinearAlgebra.Complex
 {
+    using Distributions;
     using Generic;
     using Properties;
     using Storage;
@@ -141,6 +142,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new DiagonalMatrix(DiagonalMatrixStorage<Complex>.OfArray(array));
         }
 
+        /// <summary>
+        /// Create a new diagonal matrix and initialize each diagonal value using the provided init function.
+        /// </summary>
+        public static DiagonalMatrix Create(int rows, int columns, Func<int, Complex> init)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<Complex>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new diagonal matrix with diagonal values sampled from the provided random distribution.
+        /// </summary>
+        public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<Complex>.OfInit(rows, columns,
+                i => new Complex(distribution.Sample(), distribution.Sample())));
+        }
+
         /// <summary>
         /// Create a new diagonal matrix as a copy of the given two-dimensional array.
         /// This new matrix will be independent from the provided array.

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -136,43 +136,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static SparseMatrix Create(int rows, int columns, Func<int, int, Complex> init)
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfInit(rows, columns, init));
+        }
+
         /// <summary>
         /// Create a new sparse matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Use a dense matrix instead. Scheduled for removal in v3.0.")]
+        [Obsolete("Use a dense matrix or SparseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public SparseMatrix(int rows, int columns, Complex value)
-            : this(rows, columns)
+            : this(SparseCompressedRowMatrixStorage<Complex>.OfInit(rows, columns, (i, j) => value))
         {
-            if (value.IsZero())
-            {
-                return;
-            }
-
-            var rowPointers = _storage.RowPointers;
-            var valueCount = _storage.ValueCount = rows * columns;
-            var columnIndices = _storage.ColumnIndices = new int[valueCount];
-            var values = _storage.Values = new Complex[valueCount];
-
-            for (int i = 0, j = 0; i < values.Length; i++, j++)
-            {
-                // Reset column position to "0"
-                if (j == columns)
-                {
-                    j = 0;
-                }
-
-                values[i] = value;
-                columnIndices[i] = j;
-            }
-
-            // Set proper row pointers
-            for (var i = 0; i < rowPointers.Length; i++)
-            {
-                rowPointers[i] = ((i + 1) * columns) - columns;
-            }
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@@ -138,27 +138,40 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// Create a new dense matrix as a copy of the given enumerable.
         /// The enumerable is assumed to be in column-major order (column by column).
         /// This new matrix will be independent from the enumerable.
-        /// A new memory block will be allocated for storing the vector.
+        /// A new memory block will be allocated for storing the matrix.
         /// </summary>
         public static DenseMatrix OfColumnMajor(int rows, int columns, IEnumerable<Complex32> columnMajor)
         {
             return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static DenseMatrix Create(int rows, int columns, Func<int, int, Complex32> init)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix with values sampled from the provided random distribution.
+        /// </summary>
+        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfInit(rows, columns,
+                (i, j) => new Complex32((float) distribution.Sample(), (float) distribution.Sample())));
+        }
+
         /// <summary>
         /// Create a new dense matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Scheduled for removal in v3.0.")]
+        [Obsolete("Use DenseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public DenseMatrix(int rows, int columns, Complex32 value)
-            : this(rows, columns)
+            : this(DenseColumnMajorMatrixStorage<Complex32>.OfInit(rows, columns, (i, j) => value))
         {
-            for (var i = 0; i < _values.Length; i++)
-            {
-                _values[i] = value;
-            }
         }
 
         /// <summary>
@@ -183,19 +196,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         {
         }
 
-        /// <summary>
-        /// Create a new dense matrix with values sampled from the provided random distribution.
-        /// </summary>
-        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<Complex32>(rows, columns);
-            for (var i = 0; i < storage.Data.Length; i++)
-            {
-                storage.Data[i] = new Complex32((float)distribution.Sample(), (float)distribution.Sample());
-            }
-            return new DenseMatrix(storage);
-        }
-
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@@ -30,6 +30,7 @@
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32
 {
+    using Distributions;
     using Generic;
     using Numerics;
     using Properties;
@@ -141,6 +142,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new DiagonalMatrix(DiagonalMatrixStorage<Complex32>.OfArray(array));
         }
 
+        /// <summary>
+        /// Create a new diagonal matrix and initialize each diagonal value using the provided init function.
+        /// </summary>
+        public static DiagonalMatrix Create(int rows, int columns, Func<int, Complex32> init)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<Complex32>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new diagonal matrix with diagonal values sampled from the provided random distribution.
+        /// </summary>
+        public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<Complex32>.OfInit(rows, columns,
+                i => new Complex32((float) distribution.Sample(), (float) distribution.Sample())));
+        }
+
         /// <summary>
         /// Create a new diagonal matrix as a copy of the given two-dimensional array.
         /// This new matrix will be independent from the provided array.

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -136,43 +136,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static SparseMatrix Create(int rows, int columns, Func<int, int, Complex32> init)
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfInit(rows, columns, init));
+        }
+
         /// <summary>
         /// Create a new sparse matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Use a dense matrix instead. Scheduled for removal in v3.0.")]
+        [Obsolete("Use a dense matrix or SparseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public SparseMatrix(int rows, int columns, Complex32 value)
-            : this(rows, columns)
+            : this(SparseCompressedRowMatrixStorage<Complex32>.OfInit(rows, columns, (i,j) => value))
         {
-            if (value.IsZero())
-            {
-                return;
-            }
-
-            var rowPointers = _storage.RowPointers;
-            var valueCount = _storage.ValueCount = rows * columns;
-            var columnIndices = _storage.ColumnIndices = new int[valueCount];
-            var values = _storage.Values = new Complex32[valueCount];
-
-            for (int i = 0, j = 0; i < values.Length; i++, j++)
-            {
-                // Reset column position to "0"
-                if (j == columns)
-                {
-                    j = 0;
-                }
-
-                values[i] = value;
-                columnIndices[i] = j;
-            }
-
-            // Set proper row pointers
-            for (var i = 0; i < rowPointers.Length; i++)
-            {
-                rowPointers[i] = ((i + 1) * columns) - columns;
-            }
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@@ -140,27 +140,40 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// Create a new dense matrix as a copy of the given enumerable.
         /// The enumerable is assumed to be in column-major order (column by column).
         /// This new matrix will be independent from the enumerable.
-        /// A new memory block will be allocated for storing the vector.
+        /// A new memory block will be allocated for storing the matrix.
         /// </summary>
         public static DenseMatrix OfColumnMajor(int rows, int columns, IEnumerable<double> columnMajor)
         {
             return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static DenseMatrix Create(int rows, int columns, Func<int, int, double> init)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix with values sampled from the provided random distribution.
+        /// </summary>
+        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfInit(rows, columns,
+                (i, j) => distribution.Sample()));
+        }
+
         /// <summary>
         /// Create a new dense matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Scheduled for removal in v3.0.")]
+        [Obsolete("Use DenseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public DenseMatrix(int rows, int columns, double value)
-            : this(rows, columns)
+            : this(DenseColumnMajorMatrixStorage<double>.OfInit(rows, columns, (i, j) => value))
         {
-            for (var i = 0; i < _values.Length; i++)
-            {
-                _values[i] = value;
-            }
         }
 
         /// <summary>
@@ -185,19 +198,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
         }
 
-        /// <summary>
-        /// Create a new dense matrix with values sampled from the provided random distribution.
-        /// </summary>
-        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<double>(rows, columns);
-            for (var i = 0; i < storage.Data.Length; i++)
-            {
-                storage.Data[i] = distribution.Sample();
-            }
-            return new DenseMatrix(storage);
-        }
-
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>

src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@@ -30,6 +30,7 @@
 
 namespace MathNet.Numerics.LinearAlgebra.Double
 {
+    using Distributions;
     using Generic;
     using Properties;
     using Storage;
@@ -140,6 +141,23 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new DiagonalMatrix(DiagonalMatrixStorage<double>.OfArray(array));
         }
 
+        /// <summary>
+        /// Create a new diagonal matrix and initialize each diagonal value using the provided init function.
+        /// </summary>
+        public static DiagonalMatrix Create(int rows, int columns, Func<int, double> init)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<double>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new diagonal matrix with diagonal values sampled from the provided random distribution.
+        /// </summary>
+        public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<double>.OfInit(rows, columns,
+                i => distribution.Sample()));
+        }
+
         /// <summary>
         /// Create a new diagonal matrix as a copy of the given two-dimensional array.
         /// This new matrix will be independent from the provided array.

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -135,43 +135,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static SparseMatrix Create(int rows, int columns, Func<int, int, double> init)
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfInit(rows, columns, init));
+        }
+
         /// <summary>
         /// Create a new sparse matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Use a dense matrix instead. Scheduled for removal in v3.0.")]
+        [Obsolete("Use a dense matrix or SparseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public SparseMatrix(int rows, int columns, double value)
-            : this(rows, columns)
+            : this(SparseCompressedRowMatrixStorage<double>.OfInit(rows, columns, (i, j) => value))
         {
-            if (value == 0.0)
-            {
-                return;
-            }
-
-            var rowPointers = _storage.RowPointers;
-            var valueCount = _storage.ValueCount = rows * columns;
-            var columnIndices = _storage.ColumnIndices = new int[valueCount];
-            var values = _storage.Values = new double[valueCount];
-
-            for (int i = 0, j = 0; i < values.Length; i++, j++)
-            {
-                // Reset column position to "0"
-                if (j == columns)
-                {
-                    j = 0;
-                }
-
-                values[i] = value;
-                columnIndices[i] = j;
-            }
-
-            // Set proper row pointers
-            for (var i = 0; i < rowPointers.Length; i++)
-            {
-                rowPointers[i] = ((i + 1) * columns) - columns;
-            }
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@@ -138,27 +138,40 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// Create a new dense matrix as a copy of the given enumerable.
         /// The enumerable is assumed to be in column-major order (column by column).
         /// This new matrix will be independent from the enumerable.
-        /// A new memory block will be allocated for storing the vector.
+        /// A new memory block will be allocated for storing the matrix.
         /// </summary>
         public static DenseMatrix OfColumnMajor(int rows, int columns, IEnumerable<float> columnMajor)
         {
             return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static DenseMatrix Create(int rows, int columns, Func<int, int, float> init)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix with values sampled from the provided random distribution.
+        /// </summary>
+        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfInit(rows, columns,
+                (i, j) => (float) distribution.Sample()));
+        }
+
         /// <summary>
         /// Create a new dense matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Scheduled for removal in v3.0.")]
+        [Obsolete("Use DenseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public DenseMatrix(int rows, int columns, float value)
-            : this(rows, columns)
+            : this(DenseColumnMajorMatrixStorage<float>.OfInit(rows, columns, (i, j) => value))
         {
-            for (var i = 0; i < _values.Length; i++)
-            {
-                _values[i] = value;
-            }
         }
 
         /// <summary>
@@ -183,19 +196,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
         }
 
-        /// <summary>
-        /// Create a new dense matrix with values sampled from the provided random distribution.
-        /// </summary>
-        public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
-        {
-            var storage = new DenseColumnMajorMatrixStorage<float>(rows, columns);
-            for (var i = 0; i < storage.Data.Length; i++)
-            {
-                storage.Data[i] = (float)distribution.Sample();
-            }
-            return new DenseMatrix(storage);
-        }
-
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>

src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@@ -30,6 +30,7 @@
 
 namespace MathNet.Numerics.LinearAlgebra.Single
 {
+    using Distributions;
     using Generic;
     using Properties;
     using Storage;
@@ -140,6 +141,23 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new DiagonalMatrix(DiagonalMatrixStorage<float>.OfArray(array));
         }
 
+        /// <summary>
+        /// Create a new diagonal matrix and initialize each diagonal value using the provided init function.
+        /// </summary>
+        public static DiagonalMatrix Create(int rows, int columns, Func<int, float> init)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<float>.OfInit(rows, columns, init));
+        }
+
+        /// <summary>
+        /// Create a new diagonal matrix with diagonal values sampled from the provided random distribution.
+        /// </summary>
+        public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
+        {
+            return new DiagonalMatrix(DiagonalMatrixStorage<float>.OfInit(rows, columns,
+                i => (float) distribution.Sample()));
+        }
+
         /// <summary>
         /// Create a new diagonal matrix as a copy of the given two-dimensional array.
         /// This new matrix will be independent from the provided array.

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -135,43 +135,24 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix and initialize each value using the provided init function.
+        /// </summary>
+        public static SparseMatrix Create(int rows, int columns, Func<int, int, float> init)
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfInit(rows, columns, init));
+        }
+
         /// <summary>
         /// Create a new sparse matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to the provided value.
         /// Zero-length matrices are not supported.
         /// </summary>
         /// <exception cref="ArgumentException">If the row or column count is less than one.</exception>
-        [Obsolete("Use a dense matrix instead. Scheduled for removal in v3.0.")]
+        [Obsolete("Use a dense matrix or SparseMatrix.Create instead. Scheduled for removal in v3.0.")]
         public SparseMatrix(int rows, int columns, float value)
-            : this(rows, columns)
+            : this(SparseCompressedRowMatrixStorage<float>.OfInit(rows, columns, (i, j) => value))
         {
-            if (value == 0.0)
-            {
-                return;
-            }
-
-            var rowPointers = _storage.RowPointers;
-            var valueCount = _storage.ValueCount = rows * columns;
-            var columnIndices = _storage.ColumnIndices = new int[valueCount];
-            var values = _storage.Values = new float[valueCount];
-
-            for (int i = 0, j = 0; i < values.Length; i++, j++)
-            {
-                // Reset column position to "0"
-                if (j == columns)
-                {
-                    j = 0;
-                }
-
-                values[i] = value;
-                columnIndices[i] = j;
-            }
-
-            // Set proper row pointers
-            for (var i = 0; i < rowPointers.Length; i++)
-            {
-                rowPointers[i] = ((i + 1) * columns) - columns;
-            }
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -111,11 +111,27 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         public static DenseColumnMajorMatrixStorage<T> OfArray(T[,] array)
         {
             var storage = new DenseColumnMajorMatrixStorage<T>(array.GetLength(0), array.GetLength(1));
-            for (var i = 0; i < storage.RowCount; i++)
+            int index = 0;
+            for (var j = 0; j < storage.ColumnCount; j++)
             {
-                for (var j = 0; j < storage.ColumnCount; j++)
+                for (var i = 0; i < storage.RowCount; i++)
                 {
-                    storage.Data[(j * storage.RowCount) + i] = array[i, j];
+                    storage.Data[index++] = array[i, j];
+                }
+            }
+            return storage;
+        }
+
+        public static DenseColumnMajorMatrixStorage<T> OfInit(int rows, int columns, Func<int, int, T> init)
+        {
+            var storage = new DenseColumnMajorMatrixStorage<T>(rows, columns);
+            int index = 0;
+            for (var j = 0; j < storage.ColumnCount; j++)
+            {
+                for (var i = 0; i < storage.RowCount; i++)
+                {
+
+                    storage.Data[index++] = init(i, j);
                 }
             }
             return storage;

src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@@ -199,6 +199,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static DiagonalMatrixStorage<T> OfInit(int rows, int columns, Func<int, T> init)
+        {
+            var storage = new DiagonalMatrixStorage<T>(rows, columns);
+            for (var i = 0; i < storage.Data.Length; i++)
+            {
+                storage.Data[i] = init(i);
+            }
+            return storage;
+        }
+
         // MATRIX COPY
 
         internal override void CopyToUnchecked(MatrixStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -401,6 +401,33 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static SparseCompressedRowMatrixStorage<T> OfInit(int rows, int columns, Func<int, int, T> init)
+        {
+            var storage = new SparseCompressedRowMatrixStorage<T>(rows, columns);
+            var rowPointers = storage.RowPointers;
+            var columnIndices = new List<int>();
+            var values = new List<T>();
+
+            for (int row = 0; row < rows; row++)
+            {
+                rowPointers[row] = values.Count;
+                for (int col = 0; col < columns; col++)
+                {
+                    var item = init(row, col);
+                    if (!Zero.Equals(item))
+                    {
+                        values.Add(item);
+                        columnIndices.Add(col);
+                    }
+                }
+            }
+
+            storage.ColumnIndices = columnIndices.ToArray();
+            storage.Values = values.ToArray();
+            storage.ValueCount = values.Count;
+            return storage;
+        }
+
         public static SparseCompressedRowMatrixStorage<T> OfRowMajorEnumerable(int rows, int columns, IEnumerable<T> data)
         {
             if (data == null)