diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj
index 5a7b7f4e..61279acb 100644
--- a/src/FSharp/FSharp.fsproj
+++ b/src/FSharp/FSharp.fsproj
@@ -48,6 +48,7 @@
     <Compile Include="DenseVector.fs" />
     <Compile Include="Vector.fs" />
     <Compile Include="Main.fs" />
+    <Compile Include="Matrix.fs" />
   </ItemGroup>
   <Import Project="$(MSBuildExtensionsPath)\FSharp\1.0\Microsoft.FSharp.Targets" />
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
diff --git a/src/FSharp/Main.fs b/src/FSharp/Main.fs
index 2a6bfe62..0afbd4b8 100644
--- a/src/FSharp/Main.fs
+++ b/src/FSharp/Main.fs
@@ -33,5 +33,8 @@ open MathNet.Numerics.LinearAlgebra.Double
 /// A module which implements some F# utility functions.
 module FSharp =
 
+    /// Construct a dense matrix from a list of floating point numbers.
+    let inline matrix (lst: list<list<float>>) = DenseMatrix.of_list lst :> Matrix
+
     /// Construct a dense vector from a list of floating point numbers.
     let inline vector (lst: list<float>) = DenseVector.of_list lst :> Vector
\ No newline at end of file
diff --git a/src/FSharp/Matrix.fs b/src/FSharp/Matrix.fs
new file mode 100644
index 00000000..977aafa2
--- /dev/null
+++ b/src/FSharp/Matrix.fs
@@ -0,0 +1,198 @@
+// <copyright file="Matrix.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+
+open MathNet.Numerics.LinearAlgebra.Double
+
+/// A module which implements functional matrix operations.
+module Matrix =
+    
+    /// Fold a function over all matrix elements.
+    let inline fold (f: 'a -> float -> 'a) (acc0: 'a) (A: #Matrix) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        let mutable acc = acc0
+        for i=0 to n-1 do
+            for j=0 to m-1 do
+                acc <- f acc (A.Item(i,j))
+        acc
+
+    /// Fold a matrix by applying a given function to all matrix elements.
+    let inline foldi (f: int -> int -> 'a -> float -> 'a) (acc0: 'a) (A: #Matrix) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        let mutable acc = acc0
+        for i=0 to n-1 do
+            for j=0 to m-1 do
+                acc <- f i j acc (A.Item(i,j))
+        acc
+
+    /// Create a 2D array from a matrix.
+    let inline toArray2 (A: #Matrix) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        Array2D.init n m (fun i j -> (A.Item(i,j)))
+    
+    /// Checks whether a predicate holds for all elements of a matrix.  
+    let inline forall (p: float -> bool) (A: #Matrix) =
+        let mutable b = true
+        let mutable i = 0
+        let mutable j = 0
+        while b && i < A.RowCount do
+            b <- b && (p (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Chechks whether a predicate holds for at least one element of a matrix.
+    let inline exists (p: float -> bool) (A: #Matrix) =
+        let mutable b = false
+        let mutable i = 0
+        let mutable j = 0
+        while not(b) && i < A.RowCount do
+            b <- b || (p (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Checks whether a position dependent predicate holds for all elements of a matrix.
+    let inline foralli (p: int -> int -> float -> bool) (A: #Matrix) =
+        let mutable b = true
+        let mutable i = 0
+        let mutable j = 0
+        while b && i < A.RowCount do
+            b <- b && (p i j (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Checks whether a position dependent predicate holds for at least one element of a matrix.
+    let inline existsi (p: int -> int -> float -> bool) (A: #Matrix) =
+        let mutable b = false
+        let mutable i = 0
+        let mutable j = 0
+        while not(b) && i < A.RowCount do
+            b <- b || (p i j (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Map every matrix element using the given function.
+    let inline map (f: float -> float) (A: #Matrix) =
+        let N = A.RowCount
+        let M = A.ColumnCount
+        let C = A.Clone()
+        for i=0 to N-1 do
+            for j=0 to M-1 do
+                C.[i,j] <- f (C.Item(i,j))
+        C
+    
+    /// Map every matrix element using the given position dependent function.
+    let inline mapi (f: int -> int -> float -> float) (A: #Matrix) =
+        let N = A.RowCount
+        let M = A.ColumnCount
+        let C = A.Clone()
+        for i=0 to N-1 do
+            for j=0 to M-1 do
+                C.[i,j] <- f i j (C.Item(i,j))
+        C
+    
+    /// In-place assignment.
+    let inline inplaceAssign (f: int -> int -> float) (A: #Matrix) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                A.Item(i,j) <- f i j
+    
+    /// In-place map of every matrix element using a position dependent function.
+    let inline inplaceMapi (f: int -> int -> float -> float) (A: #Matrix) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                A.Item(i,j) <- f i j (A.Item(i,j))
+        
+    /// Creates a sequence that iterates the non-zero entries in the matrix.
+    let inline nonZeroEntries (A: #Matrix) =
+        seq { for i in 0 .. A.RowCount-1 do
+                for j in 0 .. A.ColumnCount-1 do
+                  if A.Item(i,j) <> 0.0 then yield (i,j, A.Item(i,j)) }
+    
+    /// Returns the sum of all elements of a matrix.
+    let inline sum (A: #Matrix) =
+        let mutable f = 0.0
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f <- f + A.Item(i,j)
+        f
+    
+    /// Iterates over all elements of a matrix.
+    let inline iter (f: float -> unit) (A: #Matrix) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f (A.Item(i,j))
+        ()
+    
+    /// Iterates over all elements of a matrix using the element indices.
+    let inline iteri (f: int -> int -> float -> unit) (A: #Matrix) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f i j (A.Item(i,j))
+        ()
+        
+    /// Fold one column.
+    let inline foldCol (f: 'a -> float -> 'a) acc (A: #Matrix) k =
+        let mutable macc = acc
+        for i=0 to A.RowCount-1 do
+            macc <- f macc (A.Item(i,k))
+        macc
+    
+    /// Fold one row.
+    let inline foldRow (f: 'a -> float -> 'a) acc (A: #Matrix) k =
+        let mutable macc = acc
+        for i=0 to A.ColumnCount-1 do
+            macc <- f macc (A.Item(k,i))
+        macc
+
+    /// Fold all columns into one row vector.
+    let inline foldByCol (f: float -> float -> float) acc (A: #Matrix) =
+        let v = new DenseVector(A.ColumnCount)
+        for k=0 to A.ColumnCount-1 do
+            let mutable macc = acc
+            for i=0 to A.RowCount-1 do
+                macc <- f macc (A.Item(i,k))
+            v.[k] <- macc
+        v :> Vector
+    
+    /// Fold all rows into one column vector.
+    let inline foldByRow (f: float -> float -> float) acc (A: #Matrix) =
+        let v = new DenseVector(A.RowCount)
+        for k=0 to A.RowCount-1 do
+            let mutable macc = acc
+            for i=0 to A.ColumnCount-1 do
+                macc <- f macc (A.Item(k,i))
+            v.[k] <- macc
+        v :> Vector
\ No newline at end of file
diff --git a/src/FSharpUnitTests/Program.fs b/src/FSharpUnitTests/Program.fs
index 3d1e98d3..c088b4e4 100644
--- a/src/FSharpUnitTests/Program.fs
+++ b/src/FSharpUnitTests/Program.fs
@@ -1,27 +1,78 @@
 open FsUnit
 open MathNet.Numerics.FSharp
-open MathNet.Numerics.LinearAlgebra
+open MathNet.Numerics.LinearAlgebra.Double
 
 /// Unit tests for the dense vector type.
 let DenseVectorTests =
 
     /// A small uniform vector.
-    let smallv = new Double.DenseVector(5, 0.3 )
+    let smallv = new DenseVector(5, 0.3 )
 
     /// A large vector with increasingly large entries
-    let largev = new Double.DenseVector( Array.init 100 (fun i -> float i / 100.0) )
+    let largev = new DenseVector( Array.init 100 (fun i -> float i / 100.0) )
     
     specs "DenseVector" [
         spec "DenseVector.init"
-            (Double.DenseVector.init 100 (fun i -> float i / 100.0) |> should equal largev)
+            (DenseVector.init 100 (fun i -> float i / 100.0) |> should equal largev)
         spec "DenseVector.of_list"
-            (Double.DenseVector.of_list [ for i in 0 .. 99 -> float i / 100.0 ] |> should equal largev)
+            (DenseVector.of_list [ for i in 0 .. 99 -> float i / 100.0 ] |> should equal largev)
         spec "DenseVector.of_seq"
-            (Double.DenseVector.of_seq (seq { for i in 0 .. 99 -> float i / 100.0 }) |> should equal largev)
+            (DenseVector.of_seq (seq { for i in 0 .. 99 -> float i / 100.0 }) |> should equal largev)
         spec "DenseVector.rangef"
-            (Double.DenseVector.rangef 0.0 0.01 0.99 |> should equal (new Double.DenseVector( [| for i in 0 .. 99 -> 0.01 * float i |] ) ))
+            (DenseVector.rangef 0.0 0.01 0.99 |> should equal (new DenseVector( [| for i in 0 .. 99 -> 0.01 * float i |] ) ))
         spec "DenseVector.range"
-            (Double.DenseVector.range 0 99 |> should equal (new Double.DenseVector( [| for i in 0 .. 99 -> float i |] ) ))
+            (DenseVector.range 0 99 |> should equal (new DenseVector( [| for i in 0 .. 99 -> float i |] ) ))
+    ]
+
+
+/// Unit tests for the matrix type.
+let MatrixTests =
+
+    /// A small uniform vector.
+    let smallM = new DenseMatrix( Array2D.create 2 2 0.3 )
+
+    /// A large vector with increasingly large entries
+    let largeM = new DenseMatrix( Array2D.init 100 100 (fun i j -> float i * 100.0 + float j) )
+    
+    specs "Matrix" [
+        spec "Matrix.fold"
+            (Matrix.fold (fun a b -> a + b) 0.0 smallM |> should equal 1.2)
+        spec "Matrix.foldi"
+            (Matrix.foldi (fun i j acc x -> acc + x + float (i+j)) 0.0 smallM |> should equal 5.2)
+        spec "Matrix.toArray2"
+            (Matrix.toArray2 smallM |> should equal (Array2D.create 2 2 0.3))
+        spec "Matrix.forall"
+            (Matrix.forall (fun x -> x = 0.3) smallM |> should equal true)
+        spec "Matrix.exists"
+            (Matrix.exists (fun x -> x = 0.5) smallM |> should equal false)
+        spec "Matrix.foralli"
+            (Matrix.foralli (fun i j x -> x = float i * 100.0 + float j) largeM |> should equal true)
+        spec "Matrix.existsi"
+            (Matrix.existsi (fun i j x -> x = float i * 100.0 + float j) largeM |> should equal true)
+        spec "Matrix.map"
+            (Matrix.map (fun x -> 2.0 * x) smallM |> should equal (2.0 * smallM))
+        spec "Matrix.mapi"
+            (Matrix.mapi (fun i j x -> float i * 100.0 + float j + x) largeM |> should equal (2.0 * largeM))
+        spec "Matrix.inplaceAssign"
+            ( let N = smallM.Clone()
+              Matrix.inplaceAssign (fun i j -> 0.0) N
+              N |> should equal (0.0 * smallM))
+        spec "Matrix.inplaceMapi"
+            ( let N = largeM.Clone()
+              Matrix.inplaceMapi (fun i j x -> 2.0 * (float i * 100.0 + float j) + x) N
+              N |> should equal (3.0 * largeM))
+        spec "Matrix.nonZeroEntries"
+            (Seq.length (Matrix.nonZeroEntries smallM) |> should equal 4)
+        spec "Matrix.sum"
+            (Matrix.sum smallM |> should equal 1.2)
+        spec "Matrix.foldCol"
+            (Matrix.foldCol (+) 0.0 largeM 0 |> should equal 495000.0)
+        spec "Matrix.foldRow"
+            (Matrix.foldRow (+) 0.0 largeM 0 |> should equal 4950.0)
+        spec "Matrix.foldByCol"
+            (Matrix.foldByCol (+) 0.0 smallM |> should equal (DenseVector.of_list [0.6;0.6] :> Vector))
+        spec "Matrix.foldByRow"
+            (Matrix.foldByRow (+) 0.0 smallM |> should equal (DenseVector.of_list [0.6;0.6] :> Vector))
     ]
     
 /// Report on errors and success and exit.
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 7c093512..52b1c4ba 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -61,7 +61,6 @@
     <Compile Include="Distributions\IContinuousDistribution.cs" />
     <Compile Include="Distributions\IDiscreteDistribution.cs" />
     <Compile Include="Distributions\IDistribution.cs" />
-    <Compile Include="Distributions\Multivariate\VectorNormal.cs" />
     <Compile Include="Distributions\Multivariate\Dirichlet.cs" />
     <Compile Include="Distributions\Multivariate\Multinomial.cs" />
     <Compile Include="IntegralTransforms\Algorithms\DiscreteHartleyTransform.Naive.cs" />