diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj
index 485af929..8c68a2e8 100644
--- a/src/FSharp/FSharp.fsproj
+++ b/src/FSharp/FSharp.fsproj
@@ -62,7 +62,6 @@
     <Compile Include="AssemblyInfo.fs" />
     <Compile Include="Random.fs" />
     <Compile Include="Distributions.fs" />
-    <Compile Include="LinearAlgebra.fs" />
     <Compile Include="LinearAlgebra.Double.Vector.fs" />
     <Compile Include="LinearAlgebra.Double.Matrix.fs" />
     <Compile Include="LinearAlgebra.Double.fs" />
@@ -74,6 +73,8 @@
     <Compile Include="Fit.fs" />
     <Compile Include="FindRoots.fs" />
     <Compile Include="RandomVariable.fs" />
+    <Compile Include="LinearAlgebra.Vector.fs" />
+    <Compile Include="LinearAlgebra.Matrix.fs" />
   </ItemGroup>
   <ItemGroup>
     <Reference Include="FSharp.Core" />
diff --git a/src/FSharp/LinearAlgebra.Double.Matrix.fs b/src/FSharp/LinearAlgebra.Double.Matrix.fs
index 2fd25abf..ee4fdd9a 100644
--- a/src/FSharp/LinearAlgebra.Double.Matrix.fs
+++ b/src/FSharp/LinearAlgebra.Double.Matrix.fs
@@ -36,161 +36,14 @@ open MathNet.Numerics.LinearAlgebra
 [<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
 module Matrix =
 
-    /// Transform a vector into a 2D array.
-    let inline toArray2 (A: #Matrix<float>) = A.ToArray()
-
-    /// In-place map of every matrix element using a function.
-    let inline mapInPlace (f: float -> float) (A: #Matrix<float>) =
-        A.MapInplace((fun x -> f x), true)
-
-    /// In-place map of every matrix element using a position dependent function.
-    let inline mapiInPlace (f: int -> int -> float -> float) (A: #Matrix<float>) =
-        A.MapIndexedInplace((fun i j x -> f i j x), true)
-
-    /// In-place map of every matrix element using a function.
-    /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapnzInPlace (f: float -> float) (A: #Matrix<float>) =
-        A.MapInplace((fun x -> f x), false)
-
-    /// In-place map of every matrix element using a position dependent function.
-    /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapinzInPlace (f: int -> int -> float -> float) (A: #Matrix<float>) =
-        A.MapIndexedInplace((fun i j x -> f i j x), false)
-
-    /// In-place map every matrix column using the given position dependent function.
-    let inline mapColsInPlace (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        for j = 0 to A.ColumnCount-1 do
-            A.SetColumn(j, f j (A.Column(j)))
-
-    /// In-place map every matrix row using the given position dependent function.
-    let inline mapRowsInPlace (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        for i = 0 to A.RowCount-1 do
-            A.SetRow(i, f i (A.Row(i)))
-
-    /// Map every matrix element using the given function.
-    let inline map (f: float -> float) (A: #Matrix<float>) =
-        let A = A.Clone()
-        A.MapInplace((fun x -> f x), true)
-        A
-
-    /// Map every matrix element using the given function.
-    /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapnz (f: float -> float) (A: #Matrix<float>) =
-        let A = A.Clone()
-        A.MapInplace((fun x -> f x), false)
-        A
-
-    /// Map every matrix element using the given position dependent function.
-    let inline mapi (f: int -> int -> float -> float) (A: #Matrix<float>) =
-        let A = A.Clone()
-        A.MapIndexedInplace((fun i j x -> f i j x), true)
-        A
-
-    /// Map every matrix element using the given position dependent function.
-    /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapinz (f: int -> int -> float -> float) (A: #Matrix<float>) =
-        let A = A.Clone()
-        A.MapIndexedInplace((fun i j x -> f i j x), false)
-        A
-
-    /// Map every matrix column using the given position dependent function.
-    let inline mapCols (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        let A = A.Clone()
-        mapColsInPlace f A
-        A
-
-    /// Map every matrix row using the given position dependent function.
-    let inline mapRows (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        let A = A.Clone()
-        mapRowsInPlace f A
-        A
-
-    /// Fold a function over all matrix elements.
-    let inline fold (f: 'a -> float -> 'a) (acc0: 'a) (A: #Matrix<float>) =
-        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.At(i,j))
-        acc
-
-    /// Fold a function over all matrix elements in reverse order.
-    let inline foldBack (f: float -> 'a -> 'a) (acc0: 'a) (A: #Matrix<float>) =
-        let n = A.RowCount
-        let m = A.ColumnCount
-        let mutable acc = acc0
-        for i in n-1 .. -1 .. 0 do
-            for j in m-1 .. -1 .. 0 do
-                acc <- f (A.At(i,j)) acc
-        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<float>) =
-        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.At(i,j))
-        acc
-
-    /// Checks whether a predicate holds for all elements of a matrix.
-    let inline forall (p: float -> bool) (A: #Matrix<float>) =
-        let mutable b = true
-        let mutable i = 0
-        let mutable j = 0
-        while b && i < A.RowCount do
-            b <- b && (p (A.At(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<float>) =
-        let mutable b = false
-        let mutable i = 0
-        let mutable j = 0
-        while not(b) && i < A.RowCount do
-            b <- b || (p (A.At(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<float>) =
-        let mutable b = true
-        let mutable i = 0
-        let mutable j = 0
-        while b && i < A.RowCount do
-            b <- b && (p i j (A.At(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<float>) =
-        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.At(i,j)))
-            j <- j+1
-            if j = A.ColumnCount then i <- i+1; j <- 0
-        b
-
-    /// In-place assignment.
-    let inline inplaceAssign (f: int -> int -> float) (A: #Matrix<float>) =
-        A.MapIndexedInplace((fun i j x -> f i j), true)
-
     /// Creates a sequence that iterates the non-zero entries in the matrix.
-    let inline nonZeroEntries (A: #Matrix<float>) =
+    let inline nonZeroEntries (A: #Matrix<_>) =
         seq { for i in 0 .. A.RowCount-1 do
                 for j in 0 .. A.ColumnCount-1 do
                   if A.At(i,j) <> 0.0 then yield (i, j, A.At(i,j)) }
 
     /// Returns the sum of all elements of a matrix.
-    let inline sum (A: #Matrix<float>) =
+    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
@@ -198,39 +51,13 @@ module Matrix =
         f
 
     /// Returns the sum of the results generated by applying a position dependent function to each column of the matrix.
-    let inline sumColsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<float>) =
+    let inline sumColsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<_>) =
         A.ColumnEnumerator() |> Seq.map (fun (j,col) -> f j col) |> Seq.reduce (+)
 
     /// Returns the sum of the results generated by applying a position dependent function to each row of the matrix.
-    let inline sumRowsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<float>) =
+    let inline sumRowsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<_>) =
         A.RowEnumerator() |> Seq.map (fun (i,row) -> f i row)  |> Seq.reduce (+)
 
-    /// Iterates over all elements of a matrix.
-    let inline iter (f: float -> unit) (A: #Matrix<float>) =
-        for i=0 to A.RowCount-1 do
-            for j=0 to A.ColumnCount-1 do
-                f (A.At(i,j))
-
-    /// Iterates over all elements of a matrix using the element indices.
-    let inline iteri (f: int -> int -> float -> unit) (A: #Matrix<float>) =
-        for i=0 to A.RowCount-1 do
-            for j=0 to A.ColumnCount-1 do
-                f i j (A.At(i,j))
-
-    /// Fold one column.
-    let inline foldCol (f: 'a -> float -> 'a) acc (A: #Matrix<float>) 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<float>) 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<float>) =
         let v = new DenseVector(A.ColumnCount)
diff --git a/src/FSharp/LinearAlgebra.Double.Vector.fs b/src/FSharp/LinearAlgebra.Double.Vector.fs
index 4c2606ad..a6d6b7b3 100644
--- a/src/FSharp/LinearAlgebra.Double.Vector.fs
+++ b/src/FSharp/LinearAlgebra.Double.Vector.fs
@@ -36,162 +36,6 @@ open MathNet.Numerics.LinearAlgebra
 [<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
 module Vector =
 
-    /// Transform a vector into an array.
-    let inline toArray (v: #Vector<float>) = v.ToArray()
-
-    /// Transform a vector into a list.
-    let inline toList (v: #Vector<float>) = List.init v.Count v.At
-
-    /// In-place mutation by applying a function to every element of the vector.
-    let inline mapInPlace (f: float -> float) (v: #Vector<float>) =
-        v.MapInplace((fun x -> f x), true)
-
-    /// In-place mutation by applying a function to every element of the vector.
-    let inline mapiInPlace (f: int -> float -> float) (v: #Vector<float>) =
-        v.MapIndexedInplace((fun i x -> f i x), true)
-
-    /// In-place mutation by applying a function to every element of the vector.
-    /// Zero-values may be skipped (relevant mostly for sparse vectors).
-    let inline mapnzInPlace (f: float -> float) (v: #Vector<float>) =
-        v.MapInplace((fun x -> f x), false)
-
-    /// In-place mutation by applying a function to every element of the vector.
-    /// Zero-values may be skipped (relevant mostly for sparse vectors).
-    let inline mapinzInPlace (f: int -> float -> float) (v: #Vector<float>) =
-        v.MapIndexedInplace((fun i x -> f i x), false)
-
-    /// Maps a vector to a new vector by applying a function to every element.
-    let inline map f (v: #Vector<float>) =
-        let w = v.Clone()
-        w.MapInplace((fun x -> f x), true)
-        w
-
-    /// Maps a vector to a new vector by applying a function to every element.
-    /// Zero-values may be skipped (relevant mostly for sparse vectors).
-    let inline mapnz f (v: #Vector<float>) =
-        let w = v.Clone()
-        w.MapInplace((fun x -> f x), false)
-        w
-
-    /// Maps a vector to a new vector by applying a function to every element.
-    let inline mapi (f: int -> float -> float) (v: #Vector<float>) =
-        let w = v.Clone()
-        w.MapIndexedInplace((fun i x -> f i x), true)
-        w
-
-    /// Maps a vector to a new vector by applying a function to every element.
-    /// Zero-values may be skipped (relevant mostly for sparse vectors).
-    let inline mapinz (f: int -> float -> float) (v: #Vector<float>) =
-        let w = v.Clone()
-        w.MapIndexedInplace((fun i x -> f i x), false)
-        w
-
-    /// In-place vector addition.
-    let inline addInPlace (v: #Vector<float>) (w: #Vector<float>) = v.Add(w, v)
-
-    /// In place vector subtraction.
-    let inline subInPlace (v: #Vector<float>) (w: #Vector<float>) = v.Subtract(w, v)
-
-    /// Applies a function to all elements of the vector.
-    let inline iter (f: float -> unit) (v: #Vector<float>) =
-        for i=0 to v.Count-1 do
-            f (v.At i)
-
-    /// Applies a function to all elements of the vector.
-    let inline iteri (f: int -> float -> unit) (v: #Vector<float>) =
-        for i=0 to v.Count-1 do
-            f i (v.At i)
-
-
-    /// Fold all entries of a vector.
-    let inline fold (f: 'a -> float -> 'a) (acc0: 'a) (v: #Vector<float>) =
-        let mutable acc = acc0
-        for i=0 to v.Count-1 do
-            acc <- f acc (v.At i)
-        acc
-
-    /// Fold all entries of a vector in reverse order.
-    let inline foldBack (f: float -> 'a -> 'a) (acc0: 'a) (v: #Vector<float>) =
-        let mutable acc = acc0
-        for i=2 to v.Count do
-            acc <- f (v.At (v.Count - i)) acc
-        acc
-
-    /// Fold all entries of a vector using a position dependent folding function.
-    let inline foldi (f: int -> 'a -> float -> 'a) (acc0: 'a) (v: #Vector<float>) =
-        let mutable acc = acc0
-        for i=0 to v.Count-1 do
-            acc <- f i acc (v.At i)
-        acc
-
-    /// Checks whether a predicate is satisfied for every element in the vector.
-    let inline forall (p: float -> bool) (v: #Vector<float>) =
-        let mutable b = true
-        let mutable i = 0
-        while b && i < v.Count do
-            b <- b && (p (v.At i))
-            i <- i+1
-        b
-
-    /// Checks whether there is an entry in the vector that satisfies a given predicate.
-    let inline exists (p: float -> bool) (v: #Vector<float>) =
-        let mutable b = false
-        let mutable i = 0
-        while not(b) && i < v.Count do
-            b <- b || (p (v.At i))
-            i <- i+1
-        b
-
-    /// Checks whether a predicate is true for all entries in a vector.
-    let inline foralli (p: int -> float -> bool) (v: #Vector<float>) =
-        let mutable b = true
-        let mutable i = 0
-        while b && i < v.Count do
-            b <- b && (p i (v.At i))
-            i <- i+1
-        b
-
-    /// Checks whether there is an entry in the vector that satisfies a given position dependent predicate.
-    let inline existsi (p: int -> float -> bool) (v: #Vector<float>) =
-        let mutable b = false
-        let mutable i = 0
-        while not(b) && i < v.Count do
-            b <- b || (p i (v.At i))
-            i <- i+1
-        b
-
-    /// Scans a vector; like fold but returns the intermediate result.
-    let inline scan (f: float -> float -> float) (v: #Vector<float>) =
-        let w = v.Clone()
-        let mutable p = v.Item(0)
-        for i=1 to v.Count-1 do
-            p <- f p (v.At i)
-            w.At(i, p)
-        w
-
-    /// Scans a vector in reverse order; like foldBack but returns the intermediate result.
-    let inline scanBack (f: float -> float -> float) (v: #Vector<float>) =
-        let w = v.Clone()
-        let mutable p = v.At (v.Count-1)
-        for i=2 to v.Count do
-            p <- f (v.At (v.Count - i)) p
-            w.At(v.Count - i, p)
-        w
-
-    /// Reduces a vector: the result of this function will be f(...f(f(v[0],v[1]), v[2]),..., v[n]).
-    let inline reduce (f: float -> float -> float) (v: #Vector<float>) =
-        let mutable p = v.Item(0)
-        for i=1 to v.Count-1 do
-            p <- f p (v.At i)
-        p
-
-    /// Reduces a vector in reverse order: the result of this function will be f(v[1], ..., f(v[n-2], f(v[n-1],v[n]))...).
-    let inline reduceBack (f: float -> float -> float) (v: #Vector<float>) =
-        let mutable p = v.Item(v.Count-1)
-        for i=2 to v.Count do
-            p <- f (v.At (v.Count - i)) p
-        p
-
     /// Creates a new vector and inserts the given value at the given index.
     let inline insert index value (v: #Vector<float>) =
         let newV = new DenseVector(v.Count + 1)
diff --git a/src/FSharp/LinearAlgebra.Double.fs b/src/FSharp/LinearAlgebra.Double.fs
index 77dd36cd..5f3adb33 100644
--- a/src/FSharp/LinearAlgebra.Double.fs
+++ b/src/FSharp/LinearAlgebra.Double.fs
@@ -37,7 +37,7 @@ open MathNet.Numerics.LinearAlgebra
 module Utility =
 
     /// Construct a dense matrix from a list of floating point numbers.
-    let inline matrix (lst: list<list<float>>) = DenseMatrix.ofList lst :> Matrix<float>
+    let inline matrix (lst: list<list<float>>) = DenseMatrix.ofList lst
 
     /// Construct a dense vector from a list of floating point numbers.
-    let inline vector (lst: list<float>) = DenseVector.ofList lst :> Vector<float>
+    let inline vector (lst: list<float>) = DenseVector.ofList lst
diff --git a/src/FSharp/LinearAlgebra.Matrix.fs b/src/FSharp/LinearAlgebra.Matrix.fs
new file mode 100644
index 00000000..d2846fd5
--- /dev/null
+++ b/src/FSharp/LinearAlgebra.Matrix.fs
@@ -0,0 +1,268 @@
+// <copyright file="LinearAlgebra.Matrix.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 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
+
+/// Module that contains implementation of useful F#-specific extension members for generic matrices
+[<AutoOpen>]
+module MatrixExtensions =
+
+  // A type extension for the generic matrix type that
+  // adds the 'GetSlice' method to allow m.[r1 .. r2, c1 .. c2] syntax
+  type MathNet.Numerics.LinearAlgebra.
+      Matrix<'T when 'T : struct and 'T : (new : unit -> 'T)
+                 and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable
+                 and 'T :> System.ValueType> with
+
+    /// Gets a submatrix using a specified column range and
+    /// row range (all indices are optional)
+    /// This method can be used via the x.[r1 .. r2, c1 .. c2 ] syntax
+    member x.GetSlice(rstart, rfinish, cstart, cfinish) =
+      let cstart = defaultArg cstart 0
+      let rstart = defaultArg rstart 0
+      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
+      let rfinish = defaultArg rfinish (x.RowCount - 1)
+      x.SubMatrix(rstart, rfinish - rstart + 1, cstart, cfinish - cstart + 1)
+
+    /// Sets a submatrix using a specified column range and
+    /// row range (all indices are optional)
+    /// This method can be used via the x.[r1 .. r2, c1 .. c2 ] <- m syntax
+    member x.SetSlice(rstart, rfinish, cstart, cfinish, values) =
+      let cstart = defaultArg cstart 0
+      let rstart = defaultArg rstart 0
+      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
+      let rfinish = defaultArg rfinish (x.RowCount - 1)
+      x.SetSubMatrix(rstart, rfinish - rstart + 1, cstart, cfinish - cstart + 1, values)
+
+    /// Gets a row subvector using a specified row index and column range.
+    /// This method can be used via the x.[r, c1 .. c2] syntax (F#3.1)
+    member x.GetSlice(r, cstart, cfinish) =
+      let cstart = defaultArg cstart 0
+      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
+      x.Row(r, cstart, cfinish - cstart + 1)
+
+    /// Gets a column subvector using a specified row index and column range.
+    /// This method can be used via the x.[r1 .. r2, c] syntax (F#3.1)
+    member x.GetSlice(rstart, rfinish, c) =
+      let rstart = defaultArg rstart 0
+      let rfinish = defaultArg rfinish (x.RowCount - 1)
+      x.Column(c, rstart, rfinish - rstart + 1)
+
+    /// Sets a row subvector using a specified row index and column range.
+    /// This method can be used via the x.[r, c1 .. c2] <- v syntax (F#3.1)
+    member x.SetSlice(r, cstart, cfinish, values) =
+      let cstart = defaultArg cstart 0
+      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
+      x.SetRow(r, cstart, cfinish - cstart + 1, values)
+
+    /// Sets a column subvector using a specified row index and column range.
+    /// This method can be used via the x.[r1 .. r2, c] <- v syntax (F#3.1)
+    member x.SetSlice(rstart, rfinish, c, values) =
+      let rstart = defaultArg rstart 0
+      let rfinish = defaultArg rfinish (x.RowCount - 1)
+      x.SetColumn(c, rstart, rfinish - rstart + 1, values)
+
+
+/// A module which implements functional matrix operations.
+[<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
+module Matrix =
+
+    /// Transform a vector into a 2D array.
+    let inline toArray2 (A: #Matrix<_>) = A.ToArray()
+
+    /// In-place map of every matrix element using a function.
+    let inline mapInPlace f (A: #Matrix<_>) =
+        A.MapInplace((fun x -> f x), true)
+
+    /// In-place map of every matrix element using a position dependent function.
+    let inline mapiInPlace f (A: #Matrix<_>) =
+        A.MapIndexedInplace((fun i j x -> f i j x), true)
+
+    /// In-place map of every matrix element using a function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapnzInPlace f (A: #Matrix<_>) =
+        A.MapInplace((fun x -> f x), false)
+
+    /// In-place map of every matrix element using a position dependent function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapinzInPlace f (A: #Matrix<_>) =
+        A.MapIndexedInplace((fun i j x -> f i j x), false)
+
+    /// In-place map every matrix column using the given position dependent function.
+    let inline mapColsInPlace (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) =
+        for j = 0 to A.ColumnCount-1 do
+            A.SetColumn(j, f j (A.Column(j)))
+
+    /// In-place map every matrix row using the given position dependent function.
+    let inline mapRowsInPlace (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) =
+        for i = 0 to A.RowCount-1 do
+            A.SetRow(i, f i (A.Row(i)))
+
+    /// Map every matrix element using the given function.
+    let inline map f (A: #Matrix<_>) =
+        let A = A.Clone()
+        A.MapInplace((fun x -> f x), true)
+        A
+
+    /// Map every matrix element using the given function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapnz f (A: #Matrix<_>) =
+        let A = A.Clone()
+        A.MapInplace((fun x -> f x), false)
+        A
+
+    /// Map every matrix element using the given position dependent function.
+    let inline mapi f (A: #Matrix<_>) =
+        let A = A.Clone()
+        A.MapIndexedInplace((fun i j x -> f i j x), true)
+        A
+
+    /// Map every matrix element using the given position dependent function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapinz f (A: #Matrix<_>) =
+        let A = A.Clone()
+        A.MapIndexedInplace((fun i j x -> f i j x), false)
+        A
+
+    /// Map every matrix column using the given position dependent function.
+    let inline mapCols (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) =
+        let A = A.Clone()
+        mapColsInPlace f A
+        A
+
+    /// Map every matrix row using the given position dependent function.
+    let inline mapRows (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) =
+        let A = A.Clone()
+        mapRowsInPlace f A
+        A
+
+    /// Fold a function over all matrix elements.
+    let inline fold f acc0 (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.At(i,j))
+        acc
+
+    /// Fold a function over all matrix elements in reverse order.
+    let inline foldBack f acc0 (A: #Matrix<_>) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        let mutable acc = acc0
+        for i in n-1 .. -1 .. 0 do
+            for j in m-1 .. -1 .. 0 do
+                acc <- f (A.At(i,j)) acc
+        acc
+
+    /// Fold a matrix by applying a given function to all matrix elements.
+    let inline foldi f acc0 (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.At(i,j))
+        acc
+
+    /// Checks whether a predicate holds for all elements of a matrix.
+    let inline forall p (A: #Matrix<_>) =
+        let mutable b = true
+        let mutable i = 0
+        let mutable j = 0
+        while b && i < A.RowCount do
+            b <- b && (p (A.At(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 (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.At(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 (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.At(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 (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.At(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+
+    /// In-place assignment.
+    let inline inplaceAssign f (A: #Matrix<_>) =
+        A.MapIndexedInplace((fun i j x -> f i j), true)
+
+    /// Iterates over all elements of a matrix.
+    let inline iter f (A: #Matrix<_>) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f (A.At(i,j))
+
+    /// Iterates over all elements of a matrix using the element indices.
+    let inline iteri f (A: #Matrix<_>) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f i j (A.At(i,j))
+
+    /// Fold one column.
+    let inline foldCol f 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 acc (A: #Matrix<_>) k =
+        let mutable macc = acc
+        for i=0 to A.ColumnCount-1 do
+            macc <- f macc (A.Item(k,i))
+        macc
diff --git a/src/FSharp/LinearAlgebra.Vector.fs b/src/FSharp/LinearAlgebra.Vector.fs
new file mode 100644
index 00000000..1b2ebf60
--- /dev/null
+++ b/src/FSharp/LinearAlgebra.Vector.fs
@@ -0,0 +1,220 @@
+// <copyright file="LinearAlgebra.Vector.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 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
+
+/// Module that contains implementation of useful F#-specific extension members for generic vectors
+[<AutoOpen>]
+module VectorExtensions =
+
+  // A type extension for the generic vector type that
+  // adds the 'GetSlice' method to allow vec.[a .. b] syntax
+  type MathNet.Numerics.LinearAlgebra.
+      Vector<'T when 'T : struct and 'T : (new : unit -> 'T)
+                 and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable
+                 and 'T :> System.ValueType> with
+
+    /// Gets a slice of a vector starting at a specified index
+    /// and ending at a specified index (both indices are optional)
+    /// This method can be used via the x.[start .. finish] syntax
+    member x.GetSlice(start, finish) =
+      let start = defaultArg start 0
+      let finish = defaultArg finish (x.Count - 1)
+      x.SubVector(start, finish - start + 1)
+
+    /// Sets a slice of a vector starting at a specified index
+    /// and ending at a specified index (both indices are optional)
+    /// This method can be used via the x.[start .. finish] <- v syntax
+    member x.SetSlice(start, finish, values) =
+      let start = defaultArg start 0
+      let finish = defaultArg finish (x.Count - 1)
+      x.SetSubVector(start, finish - start + 1, values)
+
+
+/// A module which implements functional vector operations.
+[<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
+module Vector =
+
+    /// Transform a vector into an array.
+    let inline toArray (v: #Vector<_>) = v.ToArray()
+
+    /// Transform a vector into a list.
+    let inline toList (v: #Vector<_>) = List.init v.Count v.At
+
+    /// In-place mutation by applying a function to every element of the vector.
+    let inline mapInPlace f (v: #Vector<_>) =
+        v.MapInplace((fun x -> f x), true)
+
+    /// In-place mutation by applying a function to every element of the vector.
+    let inline mapiInPlace f (v: #Vector<_>) =
+        v.MapIndexedInplace((fun i x -> f i x), true)
+
+    /// In-place mutation by applying a function to every element of the vector.
+    /// Zero-values may be skipped (relevant mostly for sparse vectors).
+    let inline mapnzInPlace f (v: #Vector<_>) =
+        v.MapInplace((fun x -> f x), false)
+
+    /// In-place mutation by applying a function to every element of the vector.
+    /// Zero-values may be skipped (relevant mostly for sparse vectors).
+    let inline mapinzInPlace (f: int -> float -> float) (v: #Vector<float>) =
+        v.MapIndexedInplace((fun i x -> f i x), false)
+
+    /// Maps a vector to a new vector by applying a function to every element.
+    let inline map f (v: #Vector<_>) =
+        let w = v.Clone()
+        w.MapInplace((fun x -> f x), true)
+        w
+
+    /// Maps a vector to a new vector by applying a function to every element.
+    /// Zero-values may be skipped (relevant mostly for sparse vectors).
+    let inline mapnz f (v: #Vector<_>) =
+        let w = v.Clone()
+        w.MapInplace((fun x -> f x), false)
+        w
+
+    /// Maps a vector to a new vector by applying a function to every element.
+    let inline mapi f (v: #Vector<_>) =
+        let w = v.Clone()
+        w.MapIndexedInplace((fun i x -> f i x), true)
+        w
+
+    /// Maps a vector to a new vector by applying a function to every element.
+    /// Zero-values may be skipped (relevant mostly for sparse vectors).
+    let inline mapinz f (v: #Vector<_>) =
+        let w = v.Clone()
+        w.MapIndexedInplace((fun i x -> f i x), false)
+        w
+
+    /// In-place vector addition.
+    let inline addInPlace (v: #Vector<_>) (w: #Vector<_>) = v.Add(w, v)
+
+    /// In place vector subtraction.
+    let inline subInPlace (v: #Vector<_>) (w: #Vector<_>) = v.Subtract(w, v)
+
+    /// Applies a function to all elements of the vector.
+    let inline iter f (v: #Vector<_>) =
+        for i=0 to v.Count-1 do
+            f (v.At i)
+
+    /// Applies a function to all elements of the vector.
+    let inline iteri f (v: #Vector<_>) =
+        for i=0 to v.Count-1 do
+            f i (v.At i)
+
+
+    /// Fold all entries of a vector.
+    let inline fold f acc0 (v: #Vector<_>) =
+        let mutable acc = acc0
+        for i=0 to v.Count-1 do
+            acc <- f acc (v.At i)
+        acc
+
+    /// Fold all entries of a vector in reverse order.
+    let inline foldBack f acc0 (v: #Vector<_>) =
+        let mutable acc = acc0
+        for i=2 to v.Count do
+            acc <- f (v.At (v.Count - i)) acc
+        acc
+
+    /// Fold all entries of a vector using a position dependent folding function.
+    let inline foldi f acc0 (v: #Vector<_>) =
+        let mutable acc = acc0
+        for i=0 to v.Count-1 do
+            acc <- f i acc (v.At i)
+        acc
+
+    /// Checks whether a predicate is satisfied for every element in the vector.
+    let inline forall p (v: #Vector<_>) =
+        let mutable b = true
+        let mutable i = 0
+        while b && i < v.Count do
+            b <- b && (p (v.At i))
+            i <- i+1
+        b
+
+    /// Checks whether there is an entry in the vector that satisfies a given predicate.
+    let inline exists p (v: #Vector<_>) =
+        let mutable b = false
+        let mutable i = 0
+        while not(b) && i < v.Count do
+            b <- b || (p (v.At i))
+            i <- i+1
+        b
+
+    /// Checks whether a predicate is true for all entries in a vector.
+    let inline foralli p (v: #Vector<_>) =
+        let mutable b = true
+        let mutable i = 0
+        while b && i < v.Count do
+            b <- b && (p i (v.At i))
+            i <- i+1
+        b
+
+    /// Checks whether there is an entry in the vector that satisfies a given position dependent predicate.
+    let inline existsi p (v: #Vector<_>) =
+        let mutable b = false
+        let mutable i = 0
+        while not(b) && i < v.Count do
+            b <- b || (p i (v.At i))
+            i <- i+1
+        b
+
+    /// Scans a vector; like fold but returns the intermediate result.
+    let inline scan f (v: #Vector<_>) =
+        let w = v.Clone()
+        let mutable p = v.Item(0)
+        for i=1 to v.Count-1 do
+            p <- f p (v.At i)
+            w.At(i, p)
+        w
+
+    /// Scans a vector in reverse order; like foldBack but returns the intermediate result.
+    let inline scanBack f (v: #Vector<_>) =
+        let w = v.Clone()
+        let mutable p = v.At (v.Count-1)
+        for i=2 to v.Count do
+            p <- f (v.At (v.Count - i)) p
+            w.At(v.Count - i, p)
+        w
+
+    /// Reduces a vector: the result of this function will be f(...f(f(v[0],v[1]), v[2]),..., v[n]).
+    let inline reduce f (v: #Vector<_>) =
+        let mutable p = v.Item(0)
+        for i=1 to v.Count-1 do
+            p <- f p (v.At i)
+        p
+
+    /// Reduces a vector in reverse order: the result of this function will be f(v[1], ..., f(v[n-2], f(v[n-1],v[n]))...).
+    let inline reduceBack f (v: #Vector<_>) =
+        let mutable p = v.Item(v.Count-1)
+        for i=2 to v.Count do
+            p <- f (v.At (v.Count - i)) p
+        p
+
diff --git a/src/FSharp/LinearAlgebra.fs b/src/FSharp/LinearAlgebra.fs
deleted file mode 100644
index e348f464..00000000
--- a/src/FSharp/LinearAlgebra.fs
+++ /dev/null
@@ -1,115 +0,0 @@
-// <copyright file="LinearAlgebra.fs" company="Math.NET">
-// Math.NET Numerics, part of the Math.NET Project
-// http://numerics.mathdotnet.com
-// http://github.com/mathnet/mathnet-numerics
-// http://mathnetnumerics.codeplex.com
-//
-// Copyright (c) 2009-2013 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
-
-// Module that contains implementation of useful F#-specific
-// extension members for generic Matrix and Vector types
-[<AutoOpen>]
-module FSharpExtensions =
-
-  // A type extension for the generic vector type that
-  // adds the 'GetSlice' method to allow vec.[a .. b] syntax
-  type MathNet.Numerics.LinearAlgebra.
-      Vector<'T when 'T : struct and 'T : (new : unit -> 'T)
-                 and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable
-                 and 'T :> System.ValueType> with
-
-    /// Gets a slice of a vector starting at a specified index
-    /// and ending at a specified index (both indices are optional)
-    /// This method can be used via the x.[start .. finish] syntax
-    member x.GetSlice(start, finish) =
-      let start = defaultArg start 0
-      let finish = defaultArg finish (x.Count - 1)
-      x.SubVector(start, finish - start + 1)
-
-    /// Sets a slice of a vector starting at a specified index
-    /// and ending at a specified index (both indices are optional)
-    /// This method can be used via the x.[start .. finish] <- v syntax
-    member x.SetSlice(start, finish, values) =
-      let start = defaultArg start 0
-      let finish = defaultArg finish (x.Count - 1)
-      x.SetSubVector(start, finish - start + 1, values)
-
-
-  // A type extension for the generic matrix type that
-  // adds the 'GetSlice' method to allow m.[r1 .. r2, c1 .. c2] syntax
-  type MathNet.Numerics.LinearAlgebra.
-      Matrix<'T when 'T : struct and 'T : (new : unit -> 'T)
-                 and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable
-                 and 'T :> System.ValueType> with
-
-    /// Gets a submatrix using a specified column range and
-    /// row range (all indices are optional)
-    /// This method can be used via the x.[r1 .. r2, c1 .. c2 ] syntax
-    member x.GetSlice(rstart, rfinish, cstart, cfinish) =
-      let cstart = defaultArg cstart 0
-      let rstart = defaultArg rstart 0
-      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
-      let rfinish = defaultArg rfinish (x.RowCount - 1)
-      x.SubMatrix(rstart, rfinish - rstart + 1, cstart, cfinish - cstart + 1)
-
-    /// Sets a submatrix using a specified column range and
-    /// row range (all indices are optional)
-    /// This method can be used via the x.[r1 .. r2, c1 .. c2 ] <- m syntax
-    member x.SetSlice(rstart, rfinish, cstart, cfinish, values) =
-      let cstart = defaultArg cstart 0
-      let rstart = defaultArg rstart 0
-      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
-      let rfinish = defaultArg rfinish (x.RowCount - 1)
-      x.SetSubMatrix(rstart, rfinish - rstart + 1, cstart, cfinish - cstart + 1, values)
-
-    /// Gets a row subvector using a specified row index and column range.
-    /// This method can be used via the x.[r, c1 .. c2] syntax (F#3.1)
-    member x.GetSlice(r, cstart, cfinish) =
-      let cstart = defaultArg cstart 0
-      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
-      x.Row(r, cstart, cfinish - cstart + 1)
-
-    /// Gets a column subvector using a specified row index and column range.
-    /// This method can be used via the x.[r1 .. r2, c] syntax (F#3.1)
-    member x.GetSlice(rstart, rfinish, c) =
-      let rstart = defaultArg rstart 0
-      let rfinish = defaultArg rfinish (x.RowCount - 1)
-      x.Column(c, rstart, rfinish - rstart + 1)
-
-    /// Sets a row subvector using a specified row index and column range.
-    /// This method can be used via the x.[r, c1 .. c2] <- v syntax (F#3.1)
-    member x.SetSlice(r, cstart, cfinish, values) =
-      let cstart = defaultArg cstart 0
-      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
-      x.SetRow(r, cstart, cfinish - cstart + 1, values)
-
-    /// Sets a column subvector using a specified row index and column range.
-    /// This method can be used via the x.[r1 .. r2, c] <- v syntax (F#3.1)
-    member x.SetSlice(rstart, rfinish, c, values) =
-      let rstart = defaultArg rstart 0
-      let rfinish = defaultArg rfinish (x.RowCount - 1)
-      x.SetColumn(c, rstart, rfinish - rstart + 1, values)