diff --git a/src/FSharp/LinearAlgebra.Matrix.fs b/src/FSharp/LinearAlgebra.Matrix.fs
index 92a773c3..6b11f055 100644
--- a/src/FSharp/LinearAlgebra.Matrix.fs
+++ b/src/FSharp/LinearAlgebra.Matrix.fs
@@ -392,6 +392,9 @@ module DenseMatrix =
     /// Create a matrix from a 2D array of floating point numbers.
     let inline ofArray2 array = Matrix<'T>.Build.DenseOfArray(array)
 
+    /// Create a matrix from a 2D array of matrices.
+    let inline ofMatrixArray2 array = Matrix<'T>.Build.DenseOfMatrixArray(array)
+
     /// Create a matrix from a list of row vectors.
     let inline ofRows (rows: Vector<'T> list) = Matrix<'T>.Build.DenseOfRowVectors(Array.ofList rows)
 
@@ -481,6 +484,9 @@ module SparseMatrix =
     /// Create a matrix from a 2D array of floating point numbers.
     let inline ofArray2 array = Matrix<'T>.Build.SparseOfArray(array)
 
+    /// Create a matrix from a 2D array of matrices.
+    let inline ofMatrixArray2 array = Matrix<'T>.Build.SparseOfMatrixArray(array)
+
     /// Create a matrix from a list of row vectors.
     let inline ofRows (rows: Vector<'T> list) = Matrix<'T>.Build.SparseOfRowVectors(Array.ofList rows)
 
diff --git a/src/FSharpUnitTests/DenseMatrixTests.fs b/src/FSharpUnitTests/DenseMatrixTests.fs
index a791bebe..f39b59bb 100644
--- a/src/FSharpUnitTests/DenseMatrixTests.fs
+++ b/src/FSharpUnitTests/DenseMatrixTests.fs
@@ -43,6 +43,13 @@ module DenseMatrixTests =
         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.ofMatrixArray2`` () =
+        let a = DenseMatrix.ofMatrixArray2 (array2D [[smallM;smallM];[smallM;smallM];[smallM;smallM]])
+        a.[0..2,0..1] |> should equal smallM
+        a.[3..5,2..3] |> should equal smallM
+        a.[6..8,0..1] |> should equal smallM
+
     [<Test>]
     let ``DenseMatrix.ofRowSeq`` () =
         DenseMatrix.ofRowSeq (Seq.ofList [[0.3;0.3];[0.3;0.3];[0.3;0.3]]) |> should equal smallM
diff --git a/src/FSharpUnitTests/SparseMatrixTests.fs b/src/FSharpUnitTests/SparseMatrixTests.fs
index 046c2ec9..02f4c3df 100644
--- a/src/FSharpUnitTests/SparseMatrixTests.fs
+++ b/src/FSharpUnitTests/SparseMatrixTests.fs
@@ -23,6 +23,13 @@ module SparseMatrixTests =
         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.ofMatrixArray2`` () =
+        let a = SparseMatrix.ofMatrixArray2 (array2D [[smallM;smallM];[smallM;smallM];[smallM;smallM]])
+        a.[0..3,0..5] |> should equal smallM
+        a.[4..7,6..11] |> should equal smallM
+        a.[8..11,0..5] |> should equal smallM
+
     [<Test>]
     let ``SparseMatrix.ofRowSeq`` () =
         SparseMatrix.ofRowSeq (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
diff --git a/src/Numerics/LinearAlgebra/Builder.cs b/src/Numerics/LinearAlgebra/Builder.cs
index dbcf0d50..91e9431c 100644
--- a/src/Numerics/LinearAlgebra/Builder.cs
+++ b/src/Numerics/LinearAlgebra/Builder.cs
@@ -713,6 +713,39 @@ namespace MathNet.Numerics.LinearAlgebra
             return m;
         }
 
+        /// <summary>
+        /// Create a new dense matrix from a 2D array of existing matrices.
+        /// The matrices in the array are not required to be dense already.
+        /// If the matrices do not align properly, they are placed on the top left
+        /// corner of their cell with the remaining fields left zero.
+        /// </summary>
+        public Matrix<T> DenseOfMatrixArray(Matrix<T>[,] matrices)
+        {
+            var rowspans = new int[matrices.GetLength(0)];
+            var colspans = new int[matrices.GetLength(1)];
+            for (int i = 0; i < rowspans.Length; i++)
+            {
+                for (int j = 0; j < colspans.Length; j++)
+                {
+                    rowspans[i] = Math.Max(rowspans[i], matrices[i, j].RowCount);
+                    colspans[j] = Math.Max(colspans[j], matrices[i, j].ColumnCount);
+                }
+            }
+            var m = Dense(rowspans.Sum(), colspans.Sum());
+            int rowoffset = 0;
+            for (int i = 0; i < rowspans.Length; i++)
+            {
+                int coloffset = 0;
+                for (int j = 0; j < colspans.Length; j++)
+                {
+                    m.SetSubMatrix(rowoffset, coloffset, matrices[i,j]);
+                    coloffset += colspans[j];
+                }
+                rowoffset += rowspans[i];
+            }
+            return m;
+        }
+
         /// <summary>
         /// Create a new sparse matrix straight from an initialized matrix storage instance.
         /// The storage is used directly without copying.
@@ -1027,6 +1060,39 @@ namespace MathNet.Numerics.LinearAlgebra
             return m;
         }
 
+        /// <summary>
+        /// Create a new sparse matrix from a 2D array of existing matrices.
+        /// The matrices in the array are not required to be sparse already.
+        /// If the matrices do not align properly, they are placed on the top left
+        /// corner of their cell with the remaining fields left zero.
+        /// </summary>
+        public Matrix<T> SparseOfMatrixArray(Matrix<T>[,] matrices)
+        {
+            var rowspans = new int[matrices.GetLength(0)];
+            var colspans = new int[matrices.GetLength(1)];
+            for (int i = 0; i < rowspans.Length; i++)
+            {
+                for (int j = 0; j < colspans.Length; j++)
+                {
+                    rowspans[i] = Math.Max(rowspans[i], matrices[i, j].RowCount);
+                    colspans[j] = Math.Max(colspans[j], matrices[i, j].ColumnCount);
+                }
+            }
+            var m = Sparse(rowspans.Sum(), colspans.Sum());
+            int rowoffset = 0;
+            for (int i = 0; i < rowspans.Length; i++)
+            {
+                int coloffset = 0;
+                for (int j = 0; j < colspans.Length; j++)
+                {
+                    m.SetSubMatrix(rowoffset, coloffset, matrices[i, j]);
+                    coloffset += colspans[j];
+                }
+                rowoffset += rowspans[i];
+            }
+            return m;
+        }
+
         /// <summary>
         /// Create a new diagonal matrix straight from an initialized matrix storage instance.
         /// The storage is used directly without copying.