LA: matrix of nested enumerables, with F# support and tests

src/FSharp/LinearAlgebra.Double.Matrix.fs

@@ -247,52 +247,59 @@ module Matrix =
 module DenseMatrix =
 
     /// Create a matrix that directly binds to a raw storage array in column-major (column by column) format, without copying.
-    let inline raw (n: int) (m: int) (columnMajor: float[]) = DenseMatrix(n, m, columnMajor)
+    let inline raw (rows: int) (cols: int) (columnMajor: float[]) = DenseMatrix(rows, cols, columnMajor)
 
     /// Create an all-zero matrix with the given dimension.
-    let inline zeroCreate (n: int) (m: int) = DenseMatrix(n, m)
+    let inline zeroCreate (rows: int) (cols: int) = DenseMatrix(rows, cols)
 
     /// Create a random matrix with the given dimension and value distribution.
-    let inline randomCreate (n: int) (m: int) dist = DenseMatrix.CreateRandom(n, m, dist)
+    let inline randomCreate (rows: int) (cols: int) dist = DenseMatrix.CreateRandom(rows, cols, dist)
 
     /// Create a matrix with the given dimension and set all values to x.
-    let inline create (n: int) (m: int) x = DenseMatrix.Create(n, m, fun i j -> x)
+    let inline create (rows: int) (cols: int) x = DenseMatrix.Create(rows, cols, fun i j -> x)
 
     /// Initialize a matrix by calling a construction function for every element.
-    let inline init (n: int) (m: int) (f: int -> int -> float) = DenseMatrix.Create(n, m, fun i j -> f i j)
+    let inline init (rows: int) (cols: int) (f: int -> int -> float) = DenseMatrix.Create(rows, cols, fun i j -> f i j)
 
     /// Create a matrix from a 2D array of floating point numbers.
     let inline ofArray2 array = DenseMatrix.OfArray(array)
 
+    /// Create a matrix from a list of sequences. Every sequence in the master sequence specifies a row.
+    /// If the dimensions are known, consider to use ofRowSeq instead to avoid multiple enumeration.
+    let inline ofSeq (fss: #seq<#seq<float>>) =
+        let n = Seq.length fss
+        let m = Seq.length (Seq.head fss)
+        DenseMatrix.OfRows(n, m, fss)
+
     /// Create a matrix from a list of float lists. Every list in the master list specifies a row.
+    /// If the dimensions are known, consider to use ofRowList instead to avoid multiple enumeration.
     let inline ofList (fll: float list list) =
         let n = List.length fll
         let m = List.length (List.head fll)
-        let A = DenseMatrix(n,m)
-        fll |> List.iteri (fun i fl ->
-                            if (List.length fl) <> m then failwith "Each subrow must be of the same length." else
-                            List.iteri (fun j f -> A.At(i,j,f)) fl)
-        A
+        DenseMatrix.OfRows(n, m, fll)
 
     /// Create a matrix from a list of sequences. Every sequence in the master sequence specifies a row.
-    let inline ofSeq (fss: #seq<#seq<float>>) =
-        let n = Seq.length fss
-        let m = Seq.length (Seq.head fss)
-        let A = DenseMatrix(n,m)
-        fss |> Seq.iteri (fun i fs ->
-                            if (Seq.length fs) <> m then failwith "Each subrow must be of the same length." else
-                            Seq.iteri (fun j f -> A.At(i,j,f)) fs)
-        A
+    let inline ofRowSeq (rows: int) (cols: int) (fss: #seq<#seq<float>>) = DenseMatrix.OfRows(rows, cols, fss)
 
-    /// Create a matrix with a given dimension from an indexed list of row, column, value tuples.
-    let inline ofListi (n: int) (m: int) (fl: list<int * int * float>) = DenseMatrix.OfIndexed(n, m, Seq.ofList fl)
+    /// Create a matrix from a list of float lists. Every list in the master list specifies a row.
+    let inline ofRowList (rows: int) (cols: int) (fll: float list list) = DenseMatrix.OfRows(rows, cols, fll)
+
+    /// Create a matrix from a list of sequences. Every sequence in the master sequence specifies a column.
+    let inline ofColumnSeq (rows: int) (cols: int) (fss: #seq<#seq<float>>) = DenseMatrix.OfColumns(rows, cols, fss)
+
+    /// Create a matrix from a list of float lists. Every list in the master list specifies a column.
+    let inline ofColumnList (rows: int) (cols: int) (fll: float list list) = DenseMatrix.OfColumns(rows, cols, fll)
 
     /// Create a matrix with a given dimension from an indexed sequences of row, column, value tuples.
-    let inline ofSeqi (n: int) (m: int) (fs: #seq<int * int * float>) = DenseMatrix.OfIndexed(n, m, fs)
+    let inline ofSeqi (rows: int) (cols: int) (fs: #seq<int * int * float>) = DenseMatrix.OfIndexed(rows, cols, fs)
+
+    /// Create a matrix with a given dimension from an indexed list of row, column, value tuples.
+    let inline ofListi (rows: int) (cols: int) (fl: list<int * int * float>) = DenseMatrix.OfIndexed(rows, cols, Seq.ofList fl)
 
     /// Create a matrix with the given entries.
-    let inline initDense (n: int) (m: int) (es: #seq<int * int * float>) =
-        let A = new DenseMatrix(n,m)
+    [<System.ObsoleteAttribute("Use ofSeqi instead. Scheduled for removal in v3.0.")>]
+    let inline initDense (rows: int) (cols: int) (es: #seq<int * int * float>) =
+        let A = new DenseMatrix(rows,cols)
         Seq.iter (fun (i,j,f) -> A.At(i,j,f)) es
         A
 
@@ -311,15 +318,15 @@ module DenseMatrix =
         A
 
     /// Initialize a matrix by calling a construction function for every row.
-    let inline initRow (n: int) (m: int) (f: int -> #Vector<float>) =
-        let A = new DenseMatrix(n,m)
-        for i=0 to n-1 do A.SetRow(i, f i)
+    let inline initRow (rows: int) (cols: int) (f: int -> #Vector<float>) =
+        let A = new DenseMatrix(rows,cols)
+        for i=0 to rows-1 do A.SetRow(i, f i)
         A
 
     /// Initialize a matrix by calling a construction function for every column.
-    let inline initCol (n: int) (m: int) (f: int -> #Vector<float>) =
-        let A = new DenseMatrix(n,m)
-        for i=0 to m-1 do A.SetColumn(i, f i)
+    let inline initCol (rows: int) (cols: int) (f: int -> #Vector<float>) =
+        let A = new DenseMatrix(rows,cols)
+        for i=0 to cols-1 do A.SetColumn(i, f i)
         A
 
 /// A module which implements functional sparse vector operations.
@@ -327,31 +334,39 @@ module DenseMatrix =
 module SparseMatrix =
 
     /// Create an all-zero matrix with the given dimension.
-    let inline zeroCreate (n: int) (m: int) = SparseMatrix(n, m)
+    let inline zeroCreate (rows: int) (cols: int) = SparseMatrix(rows, cols)
 
     /// Initialize a matrix by calling a construction function for every element.
-    let inline init (n: int) (m: int) (f: int -> int -> float) = SparseMatrix.Create(n, m, fun n m -> f n m)
+    let inline init (rows: int) (cols: int) (f: int -> int -> float) = SparseMatrix.Create(rows, cols, fun n m -> f n m)
 
     /// Create a matrix from a 2D array of floating point numbers.
     let inline ofArray2 array = SparseMatrix.OfArray(array)
 
-    /// Create a matrix from a list of float lists. Every list in the master list specifies a row.
-    let inline ofList (rows: int) (cols: int) (fll: list<int * int * float>) =
-        let A = new SparseMatrix(rows, cols)
-        fll |> List.iter (fun (i, j, x) -> A.At(i,j,x))
-        A
+    /// Create a matrix with a given dimension from an indexed sequences of row, column, value tuples.
+    [<System.ObsoleteAttribute("Use ofSeqi instead. Will be changed to expect a non-indexed seq in a future version.")>]
+    let inline ofSeq (rows: int) (cols: int) (fs: #seq<int * int * float>) = SparseMatrix.OfIndexed(rows, cols, fs)
+
+    /// Create a matrix with a given dimension from an indexed list of row, column, value tuples.
+    [<System.ObsoleteAttribute("Use ofListi instead. Will be changed to expect a non-indexed list in a future version.")>]
+    let inline ofList (rows: int) (cols: int) (fl: list<int * int * float>) = SparseMatrix.OfIndexed(rows, cols, Seq.ofList fl)
 
     /// Create a matrix from a list of sequences. Every sequence in the master sequence specifies a row.
-    let inline ofSeq (rows: int) (cols: int) (fss: #seq<int * int * float>) =
-        let A = new SparseMatrix(rows, cols)
-        fss |> Seq.iter (fun (i, j, x) -> A.At(i,j,x))
-        A
+    let inline ofRowSeq (rows: int) (cols: int) (fss: #seq<#seq<float>>) = SparseMatrix.OfRows(rows, cols, fss)
 
-    /// Create a matrix with a given dimension from an indexed list of row, column, value tuples.
-    let inline ofListi (n: int) (m: int) (fl: list<int * int * float>) = SparseMatrix.OfIndexed(n, m, Seq.ofList fl)
+    /// Create a matrix from a list of float lists. Every list in the master list specifies a row.
+    let inline ofRowList (rows: int) (cols: int) (fll: float list list) = SparseMatrix.OfRows(rows, cols, fll)
+
+    /// Create a matrix from a list of sequences. Every sequence in the master sequence specifies a column.
+    let inline ofColumnSeq (rows: int) (cols: int) (fss: #seq<#seq<float>>) = SparseMatrix.OfColumns(rows, cols, fss)
+
+    /// Create a matrix from a list of float lists. Every list in the master list specifies a column.
+    let inline ofColumnList (rows: int) (cols: int) (fll: float list list) = SparseMatrix.OfColumns(rows, cols, fll)
 
     /// Create a matrix with a given dimension from an indexed sequences of row, column, value tuples.
-    let inline ofSeqi (n: int) (m: int) (fs: #seq<int * int * float>) = SparseMatrix.OfIndexed(n, m, fs)
+    let inline ofSeqi (rows: int) (cols: int) (fs: #seq<int * int * float>) = SparseMatrix.OfIndexed(rows, cols, fs)
+
+    /// Create a matrix with a given dimension from an indexed list of row, column, value tuples.
+    let inline ofListi (rows: int) (cols: int) (fl: list<int * int * float>) = SparseMatrix.OfIndexed(rows, cols, Seq.ofList fl)
 
     /// Create a square matrix with constant diagonal entries.
     let inline constDiag (n: int) (f: float) =
@@ -368,13 +383,13 @@ module SparseMatrix =
         A
 
     /// Initialize a matrix by calling a construction function for every row.
-    let inline initRow (n: int) (m: int) (f: int -> #Vector<float>) =
-        let A = new SparseMatrix(n,m)
-        for i=0 to n-1 do A.SetRow(i, f i)
+    let inline initRow (rows: int) (cols: int) (f: int -> #Vector<float>) =
+        let A = new SparseMatrix(rows,cols)
+        for i=0 to rows-1 do A.SetRow(i, f i)
         A
 
     /// Initialize a matrix by calling a construction function for every column.
-    let inline initCol (n: int) (m: int) (f: int -> #Vector<float>) =
-        let A = new SparseMatrix(n,m)
-        for i=0 to m-1 do A.SetColumn(i, f i)
+    let inline initCol (rows: int) (cols: int) (f: int -> #Vector<float>) =
+        let A = new SparseMatrix(rows,cols)
+        for i=0 to cols-1 do A.SetColumn(i, f i)
         A

src/FSharp/LinearAlgebra.Double.Vector.fs

@@ -236,7 +236,7 @@ module SparseVector =
     let inline init (n: int) (f: int -> float) = SparseVector.Create(n, fun i -> f i)
 
     /// Create a sparse vector with a given dimension from a list of index, value pairs.
-    [<System.ObsoleteAttribute("Use ofListi instead. Will be changed to expect a non-indexed seq in a future version.")>]
+    [<System.ObsoleteAttribute("Use ofListi instead. Will be changed to expect a non-indexed list in a future version.")>]
     let inline ofList (n: int) (fl: list<int * float>) = SparseVector.OfIndexedEnumerable(n, Seq.ofList fl)
 
     /// Create a sparse vector with a given dimension from a sequence of index, value pairs.

src/FSharpUnitTests/DenseMatrixTests.fs

@@ -4,36 +4,80 @@ open NUnit.Framework
 open FsUnit
 open MathNet.Numerics.LinearAlgebra.Generic
 open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.Distributions
+open MathNet.Numerics.Statistics
 
 /// Unit tests for the dense matrix type.
 module DenseMatrixTests =
 
-    /// A small uniform vector.
-    let smallM = DenseMatrix.OfArray( Array2D.create 2 2 0.3 )
+    /// A small uniform matrix.
+    let smallM = DenseMatrix.raw 3 2 [|0.3;0.3;0.3;0.3;0.3;0.3|]
 
-    /// A large vector with increasingly large entries
-    let largeM = DenseMatrix.OfArray( Array2D.init 100 100 (fun i j -> float i * 100.0 + float j) )
+    /// A large matrix with increasingly large entries
+    let largeM =
+        Array.init (100*120) (fun k -> let (j,i) = System.Math.DivRem(k,100) in float i * 100.0 + float j)
+        |> DenseMatrix.raw 100 120
+
+    [<Test>]
+    let ``DenseMatrix.zeroCreate`` () =
+        (DenseMatrix.zeroCreate 100 120) + largeM |> should equal largeM
+
+    [<Test>]
+    let ``DenseMatrix.randomCreate`` () =
+        let m = DenseMatrix.randomCreate 100 120 (Normal.WithMeanStdDev(100.0,0.1))
+        m.Values |> ArrayStatistics.Mean |> should (equalWithin 10.0) 100.0
+        m.RowCount |> should equal 100
+        m.ColumnCount |> should equal 120
+
+    [<Test>]
+    let ``DenseMatrix.create`` () =
+        DenseMatrix.create 3 2 0.3 |> should equal smallM
 
     [<Test>]
     let ``DenseMatrix.init`` () =
-        DenseMatrix.init 100 100 (fun i j -> float i * 100.0 + float j) |> should equal largeM
+        DenseMatrix.init 3 2 (fun i j -> 0.3) |> should equal smallM
+        DenseMatrix.init 100 120 (fun i j -> float i * 100.0 + float j) |> should equal largeM
 
     [<Test>]
-    let ``DenseMatrix.ofList`` () =
-        DenseMatrix.ofList [[0.3;0.3];[0.3;0.3]] |> should equal smallM
+    let ``DenseMatrix.ofArray2`` () =
+        DenseMatrix.ofArray2 (array2D [[0.3;0.3];[0.3;0.3];[0.3;0.3]]) |> should equal smallM
+        DenseMatrix.ofArray2 (Array2D.create 3 2 0.3) |> should equal smallM
+        DenseMatrix.ofArray2 (Array2D.init 100 120 (fun i j -> float i * 100.0 + float j)) |> should equal largeM
 
     [<Test>]
     let ``DenseMatrix.ofSeq`` () =
-        DenseMatrix.ofSeq (Seq.ofList [[0.3;0.3];[0.3;0.3]]) |> should equal smallM
+        DenseMatrix.ofSeq (Seq.ofList [[0.3;0.3];[0.3;0.3];[0.3;0.3]]) |> should equal smallM
 
     [<Test>]
-    let ``DenseMatrix.ofArray2`` () =
-        DenseMatrix.ofArray2 (Array2D.create 2 2 0.3) |> should equal smallM
+    let ``DenseMatrix.ofList`` () =
+        DenseMatrix.ofList [[0.3;0.3];[0.3;0.3];[0.3;0.3]] |> should equal smallM
 
     [<Test>]
-    let ``DenseMatrix.initDense`` () =
-        DenseMatrix.initDense 100 100 (seq { for i in 0 .. 99 do
-                                             for j in 0 .. 99 -> (i,j, float i * 100.0 + float j)}) |> should equal largeM
+    let ``DenseMatrix.ofRowSeq`` () =
+        DenseMatrix.ofRowSeq 3 2 (Seq.ofList [[0.3;0.3];[0.3;0.3];[0.3;0.3]]) |> should equal smallM
+
+    [<Test>]
+    let ``DenseMatrix.ofRowList`` () =
+        DenseMatrix.ofRowList 3 2 [[0.3;0.3];[0.3;0.3];[0.3;0.3]] |> should equal smallM
+
+    [<Test>]
+    let ``DenseMatrix.ofColumnSeq`` () =
+        DenseMatrix.ofColumnSeq 3 2 (Seq.ofList [[0.3;0.3;0.3];[0.3;0.3;0.3]]) |> should equal smallM
+
+    [<Test>]
+    let ``DenseMatrix.ofColumnList`` () =
+        DenseMatrix.ofColumnList 3 2 [[0.3;0.3;0.3];[0.3;0.3;0.3]] |> should equal smallM
+
+    [<Test>]
+    let ``DenseMatrix.ofSeqi`` () =
+        seq { for i in 0 .. 99 do for j in 0 .. 119 -> (i,j, float i * 100.0 + float j)}
+        |> DenseMatrix.ofSeqi 100 120 |> should equal largeM
+
+    [<Test>]
+    let ``DenseMatrix.ofListi`` () =
+        [ for i in 0 .. 99 do for j in 0 .. 119 -> (i,j, float i * 100.0 + float j) ]
+        |> DenseMatrix.ofListi 100 120 |> should equal largeM
+
     [<Test>]
     let ``DenseMatrix.constDiag`` () =
         DenseMatrix.constDiag 100 2.0 |> should equal (2.0 * (DenseMatrix.Identity 100))
@@ -44,8 +88,8 @@ module DenseMatrixTests =
 
     [<Test>]
     let ``DenseMatrix.init_row`` () =
-        DenseMatrix.initRow 100 100 (fun i -> (DenseVector.init 100 (fun j -> float i * 100.0 + float j))) |> should equal largeM
+        DenseMatrix.initRow 100 120 (fun i -> (DenseVector.init 120 (fun j -> float i * 100.0 + float j))) |> should equal largeM
 
     [<Test>]
     let ``DenseMatrix.init_col`` () =
-        DenseMatrix.initCol 100 100 (fun j -> (DenseVector.init 100 (fun i -> float i * 100.0 + float j))) |> should equal largeM
+        DenseMatrix.initCol 100 120 (fun j -> (DenseVector.init 100 (fun i -> float i * 100.0 + float j))) |> should equal largeM

src/FSharpUnitTests/SparseMatrixTests.fs

@@ -8,16 +8,45 @@ open MathNet.Numerics.LinearAlgebra.Double
 /// Unit tests for the sparse matrix type.
 module SparseMatrixTests =
 
-    /// A small uniform vector.
-    let smallM = DenseMatrix.init 4 4 (fun i j -> if i = 1 && j = 2 then 1.0 else 0.0) :> Matrix<float>
+    /// A small uniform matrix.
+    let smallM = DenseMatrix.init 4 6 (fun i j -> if i = 1 && j = 2 then 1.0 else 0.0) :> Matrix<float>
 
     [<Test>]
-    let ``SparseMatrix.ofList`` () =
-        (SparseMatrix.ofList 4 4 [(1,2,1.0)] :> Matrix<float>) |> should equal smallM
+    let ``SparseMatrix.zeroCreate`` () =
+        (SparseMatrix.zeroCreate 4 6) + smallM |> should equal smallM
 
     [<Test>]
-    let ``SparseMatrix.ofSeq`` () =
-        (SparseMatrix.ofSeq 4 4 (Seq.ofList [(1,2,1.0)]) :> Matrix<float>) |> should equal smallM
+    let ``SparseMatrix.init`` () =
+        SparseMatrix.init 4 6 (fun i j -> if i = 1 && j = 2 then 1.0 else 0.0) |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.ofArray2`` () =
+        SparseMatrix.ofArray2 (array2D [[0.;0.;0.;0.;0.;0.];[0.;0.;1.;0.;0.;0.];[0.;0.;0.;0.;0.;0.];[0.;0.;0.;0.;0.;0.]]) |> should equal smallM
+        SparseMatrix.ofArray2 (Array2D.init 4 6 (fun i j -> if i = 1 && j = 2 then 1.0 else 0.0)) |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.ofRowSeq`` () =
+        SparseMatrix.ofRowSeq 4 6 (Seq.ofList [[0.;0.;0.;0.;0.;0.];[0.;0.;1.;0.;0.;0.];[0.;0.;0.;0.;0.;0.];[0.;0.;0.;0.;0.;0.]]) |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.ofRowList`` () =
+        SparseMatrix.ofRowList 4 6 [[0.;0.;0.;0.;0.;0.];[0.;0.;1.;0.;0.;0.];[0.;0.;0.;0.;0.;0.];[0.;0.;0.;0.;0.;0.]] |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.ofColumnSeq`` () =
+        SparseMatrix.ofColumnSeq 4 6 (Seq.ofList [[0.;0.;0.;0.];[0.;0.;0.;0.];[0.;1.;0.;0.];[0.;0.;0.;0.];[0.;0.;0.;0.];[0.;0.;0.;0.]]) |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.ofColumnList`` () =
+        SparseMatrix.ofColumnList 4 6 [[0.;0.;0.;0.];[0.;0.;0.;0.];[0.;1.;0.;0.];[0.;0.;0.;0.];[0.;0.;0.;0.];[0.;0.;0.;0.]] |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.ofSeqi`` () =
+        SparseMatrix.ofSeqi 4 6 (Seq.ofList [(1,2,1.0)]) |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.ofListi`` () =
+        SparseMatrix.ofListi 4 6 [(1,2,1.0)] |> should equal smallM
 
     [<Test>]
     let ``SparseMatrix.constDiag`` () =
@@ -26,3 +55,11 @@ module SparseMatrixTests =
     [<Test>]
     let ``SparseMatrix.diag`` () =
         SparseMatrix.diag (DenseVector.Create(100, fun i -> 2.0)) |> should equal (2.0 * (SparseMatrix.Identity 100))
+
+    [<Test>]
+    let ``SparseMatrix.init_row`` () =
+        SparseMatrix.initRow 4 6 (fun i -> if i=1 then DenseVector([|0.;0.;1.;0.;0.;0.|]) else DenseVector.zeroCreate 6) |> should equal smallM
+
+    [<Test>]
+    let ``SparseMatrix.init_col`` () =
+        SparseMatrix.initCol 4 6 (fun j -> if j=2 then DenseVector([|0.;1.;0.;0.|]) else DenseVector.zeroCreate 4) |> should equal smallM

src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@@ -156,6 +156,32 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<Complex>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<Complex>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new dense matrix and initialize each value using the provided init function.
         /// </summary>
@@ -255,7 +281,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
         /// <summary>
         /// Returns the transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The transpose of this matrix.</returns>
         public override Matrix<Complex> Transpose()
         {
@@ -287,7 +313,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         }
 
         /// <summary>Calculates the infinity norm of this matrix.</summary>
-        /// <returns>The infinity norm of this matrix.</returns>  
+        /// <returns>The infinity norm of this matrix.</returns>
         public override Complex InfinityNorm()
         {
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
@@ -589,7 +615,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The conjugate transpose of this matrix.</returns>
         public override Matrix<Complex> ConjugateTranspose()
         {
@@ -659,7 +685,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         }
 
         /// <summary>
-        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>. 
+        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>.
         /// </summary>
         /// <param name="rightSide">The matrix to get the values from.</param>
         /// <returns>A matrix containing a the same values as <paramref name="rightSide"/>.</returns>

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -147,6 +147,32 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<Complex>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<Complex>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new sparse matrix and initialize each value using the provided init function.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@@ -156,6 +156,32 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<Complex32>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<Complex32>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new dense matrix and initialize each value using the provided init function.
         /// </summary>
@@ -255,7 +281,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
         /// <summary>
         /// Returns the transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The transpose of this matrix.</returns>
         public override Matrix<Complex32> Transpose()
         {
@@ -287,7 +313,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         }
 
         /// <summary>Calculates the infinity norm of this matrix.</summary>
-        /// <returns>The infinity norm of this matrix.</returns>  
+        /// <returns>The infinity norm of this matrix.</returns>
         public override Complex32 InfinityNorm()
         {
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
@@ -589,7 +615,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The conjugate transpose of this matrix.</returns>
         public override Matrix<Complex32> ConjugateTranspose()
         {
@@ -605,7 +631,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
             return ret;
         }
-        
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
@@ -659,7 +685,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         }
 
         /// <summary>
-        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>. 
+        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>.
         /// </summary>
         /// <param name="rightSide">The matrix to get the values from.</param>
         /// <returns>A matrix containing a the same values as <paramref name="rightSide"/>.</returns>

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -147,6 +147,32 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<Complex32>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<Complex32>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new sparse matrix and initialize each value using the provided init function.
         /// </summary>

src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@@ -158,6 +158,32 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<double>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<double>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new dense matrix and initialize each value using the provided init function.
         /// </summary>
@@ -257,7 +283,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
         /// <summary>
         /// Returns the transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The transpose of this matrix.</returns>
         public override Matrix<double> Transpose()
         {
@@ -289,7 +315,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         /// <summary>Calculates the infinity norm of this matrix.</summary>
-        /// <returns>The infinity norm of this matrix.</returns>  
+        /// <returns>The infinity norm of this matrix.</returns>
         public override double InfinityNorm()
         {
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
@@ -357,7 +383,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 Control.LinearAlgebraProvider.SubtractArrays(_values, denseOther._values, denseResult._values);
             }
         }
-    
+
         /// <summary>
         /// Multiplies each element of the matrix by a scalar and places results into the result matrix.
         /// </summary>
@@ -375,7 +401,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 Control.LinearAlgebraProvider.ScaleArray(scalar, _values, denseResult._values);
             }
         }
-     
+
         /// <summary>
         /// Multiplies this matrix with a vector and places the results into the result vector.
         /// </summary>
@@ -671,7 +697,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         /// <summary>
-        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>. 
+        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>.
         /// </summary>
         /// <param name="rightSide">The matrix to get the values from.</param>
         /// <returns>A matrix containing a the same values as <paramref name="rightSide"/>.</returns>

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -146,6 +146,32 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<double>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<double>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new sparse matrix and initialize each value using the provided init function.
         /// </summary>

src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@@ -156,6 +156,32 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfColumnMajorEnumerable(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<float>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new dense matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static DenseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<float>
+        {
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new dense matrix and initialize each value using the provided init function.
         /// </summary>
@@ -255,7 +281,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
         /// <summary>
         /// Returns the transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The transpose of this matrix.</returns>
         public override Matrix<float> Transpose()
         {
@@ -287,7 +313,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         }
 
         /// <summary>Calculates the infinity norm of this matrix.</summary>
-        /// <returns>The infinity norm of this matrix.</returns>  
+        /// <returns>The infinity norm of this matrix.</returns>
         public override float InfinityNorm()
         {
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
@@ -669,7 +695,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         }
 
         /// <summary>
-        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>. 
+        /// Returns a <strong>Matrix</strong> containing the same values of <paramref name="rightSide"/>.
         /// </summary>
         /// <param name="rightSide">The matrix to get the values from.</param>
         /// <returns>A matrix containing a the same values as <paramref name="rightSide"/>.</returns>

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -68,7 +68,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
             _storage = storage;
         }
-        
+
         /// <summary>
         /// Create a new square sparse matrix with the given number of rows and columns.
         /// All cells of the matrix will be initialized to zero.
@@ -146,6 +146,32 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfColumnMajorList(rows, columns, columnMajor));
         }
 
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable columns.
+        /// Each enumerable in the master enumerable specifies a column.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfColumns<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<float>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfColumnEnumerables(rows, columns, data));
+        }
+
+        /// <summary>
+        /// Create a new sparse matrix as a copy of the given enumerable of enumerable rows.
+        /// Each enumerable in the master enumerable specifies a row.
+        /// This new matrix will be independent from the enumerables.
+        /// A new memory block will be allocated for storing the matrix.
+        /// </summary>
+        public static SparseMatrix OfRows<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<float>
+        {
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfRowEnumerables(rows, columns, data));
+        }
+
         /// <summary>
         /// Create a new sparse matrix and initialize each value using the provided init function.
         /// </summary>

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -170,6 +170,67 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return new DenseColumnMajorMatrixStorage<T>(rows, columns, array);
         }
 
+        public static DenseColumnMajorMatrixStorage<T> OfColumnEnumerables<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<T>
+        {
+            if (data == null) throw new ArgumentNullException("data");
+            var array = new T[rows*columns];
+            using (var columnIterator = data.GetEnumerator())
+            {
+                for (int column = 0; column < columns; column++)
+                {
+                    if (!columnIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
+                    var arrayColumn = columnIterator.Current as T[];
+                    if (arrayColumn != null)
+                    {
+                        Array.Copy(arrayColumn, 0, array, column*rows, rows);
+                    }
+                    else
+                    {
+                        using (var rowIterator = columnIterator.Current.GetEnumerator())
+                        {
+                            var end = (column + 1)*rows;
+                            for (int index = column*rows; index < end; index++)
+                            {
+                                if (!rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, columns));
+                                array[index] = rowIterator.Current;
+                            }
+                            if (rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, columns));
+                        }
+                    }
+                }
+                if (columnIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
+            }
+            return new DenseColumnMajorMatrixStorage<T>(rows, columns, array);
+        }
+
+        public static DenseColumnMajorMatrixStorage<T> OfRowEnumerables<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<T>
+        {
+            if (data == null) throw new ArgumentNullException("data");
+            var array = new T[rows*columns];
+            using (var rowIterator = data.GetEnumerator())
+            {
+                for (int row = 0; row < rows; row++)
+                {
+                    if (!rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
+                    using (var columnIterator = rowIterator.Current.GetEnumerator())
+                    {
+                        for (int index = row; index < array.Length; index += rows)
+                        {
+                            if (!columnIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, columns));
+                            array[index] = columnIterator.Current;
+                        }
+                        if (columnIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, columns));
+                    }
+                }
+                if (rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
+            }
+            return new DenseColumnMajorMatrixStorage<T>(rows, columns, array);
+        }
+
         // MATRIX COPY
 
         internal override void CopyToUnchecked(MatrixStorage<T> target, bool skipClearing = false)

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -471,6 +471,98 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static SparseCompressedRowMatrixStorage<T> OfRowEnumerables<TRow>(int rows, int columns, IEnumerable<TRow> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TRow : IEnumerable<T>
+        {
+            if (data == null) throw new ArgumentNullException("data");
+
+            var storage = new SparseCompressedRowMatrixStorage<T>(rows, columns);
+            var rowPointers = storage.RowPointers;
+            var columnIndices = new List<int>();
+            var values = new List<T>();
+
+            using (var rowIterator = data.GetEnumerator())
+            {
+                for (int row = 0; row < rows; row++)
+                {
+                    if (!rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
+                    rowPointers[row] = values.Count;
+                    using (var columnIterator = rowIterator.Current.GetEnumerator())
+                    {
+                        for (int col = 0; col < columns; col++)
+                        {
+                            if (!columnIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, columns));
+                            if (!Zero.Equals(columnIterator.Current))
+                            {
+                                values.Add(columnIterator.Current);
+                                columnIndices.Add(col);
+                            }
+                        }
+                        if (columnIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, columns));
+                    }
+                }
+                if (rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
+            }
+            storage.ColumnIndices = columnIndices.ToArray();
+            storage.Values = values.ToArray();
+            storage.ValueCount = values.Count;
+            return storage;
+        }
+
+        public static SparseCompressedRowMatrixStorage<T> OfColumnEnumerables<TColumn>(int rows, int columns, IEnumerable<TColumn> data)
+            // NOTE: flexible typing to 'backport' generic covariance.
+            where TColumn : IEnumerable<T>
+        {
+            if (data == null) throw new ArgumentNullException("data");
+
+            var trows = new List<Tuple<int, T>>[rows];
+            using (var columnIterator = data.GetEnumerator())
+            {
+                for (int column = 0; column < columns; column++)
+                {
+                    if (!columnIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, columns));
+                    using (var rowIterator = columnIterator.Current.GetEnumerator())
+                    {
+                        for (int row = 0; row < rows; row++)
+                        {
+                            if (!rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
+                            if (!Zero.Equals(rowIterator.Current))
+                            {
+                                var trow = trows[row] ?? (trows[row] = new List<Tuple<int, T>>());
+                                trow.Add(new Tuple<int, T>(column, rowIterator.Current));
+                            }
+                        }
+                    }
+                }
+            }
+
+            var storage = new SparseCompressedRowMatrixStorage<T>(rows, columns);
+            var rowPointers = storage.RowPointers;
+            var columnIndices = new List<int>();
+            var values = new List<T>();
+
+            int index = 0;
+            for (int row = 0; row < rows; row++)
+            {
+                rowPointers[row] = index;
+                if (trows[row] != null)
+                {
+                    foreach (var item in trows[row])
+                    {
+                        values.Add(item.Item2);
+                        columnIndices.Add(item.Item1);
+                        index++;
+                    }
+                }
+            }
+
+            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)
@@ -483,17 +575,19 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             var columnIndices = new List<int>();
             var values = new List<T>();
 
-            var iterator = data.GetEnumerator();
-            for (int row = 0; row < rows; row++)
+            using (var iterator = data.GetEnumerator())
             {
-                rowPointers[row] = values.Count;
-                for (int col = 0; col < columns; col++)
+                for (int row = 0; row < rows; row++)
                 {
-                    iterator.MoveNext();
-                    if (!Zero.Equals(iterator.Current))
+                    rowPointers[row] = values.Count;
+                    for (int col = 0; col < columns; col++)
                     {
-                        values.Add(iterator.Current);
-                        columnIndices.Add(col);
+                        iterator.MoveNext();
+                        if (!Zero.Equals(iterator.Current))
+                        {
+                            values.Add(iterator.Current);
+                            columnIndices.Add(col);
+                        }
                     }
                 }
             }