From 46ef39b84ba3ca05dabfd868bb659484e5e5b91d Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Sun, 7 Apr 2013 22:04:48 +0200
Subject: [PATCH] LA: Functional matrix init

---
 .../LinearAlgebra/Complex/DenseMatrix.cs      | 40 +++++++++----------
 .../LinearAlgebra/Complex/DiagonalMatrix.cs   | 18 +++++++++
 .../LinearAlgebra/Complex/SparseMatrix.cs     | 39 +++++-------------
 .../LinearAlgebra/Complex32/DenseMatrix.cs    | 40 +++++++++----------
 .../LinearAlgebra/Complex32/DiagonalMatrix.cs | 18 +++++++++
 .../LinearAlgebra/Complex32/SparseMatrix.cs   | 39 +++++-------------
 .../LinearAlgebra/Double/DenseMatrix.cs       | 40 +++++++++----------
 .../LinearAlgebra/Double/DiagonalMatrix.cs    | 18 +++++++++
 .../LinearAlgebra/Double/SparseMatrix.cs      | 39 +++++-------------
 .../LinearAlgebra/Single/DenseMatrix.cs       | 40 +++++++++----------
 .../LinearAlgebra/Single/DiagonalMatrix.cs    | 18 +++++++++
 .../LinearAlgebra/Single/SparseMatrix.cs      | 39 +++++-------------
 .../Storage/DenseColumnMajorMatrixStorage.cs  | 22 ++++++++--
 .../Storage/DiagonalMatrixStorage.cs          | 10 +++++
 .../SparseCompressedRowMatrixStorage.cs       | 27 +++++++++++++
 15 files changed, 248 insertions(+), 199 deletions(-)

diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index d09ded08..28d36aa8 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index 5bbbd3f2..71e0a335 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/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.
diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
index 1a99d6ca..1a96c2fe 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 74a7fec9..aa5e147a 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index ecf9d32d..6348a231 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/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.
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
index e74277c3..122e1a50 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index dd7e7590..ee262973 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 312dc04c..e812322b 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/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.
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index 0ed4c50c..ca0e3792 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 22dcb7ed..9a16110a 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index c5e6ad5e..46ea9455 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/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.
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index 6fdba1f4..f075724c 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/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>
diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
index 74e0784e..10ab5982 100644
--- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
+++ b/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;
diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
index 4980fc15..70a293f5 100644
--- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
+++ b/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)
diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
index 8dd33fe5..e13710c1 100644
--- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
+++ b/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)