diff --git a/src/FSharp/Matrix.fs b/src/FSharp/Matrix.fs
index e783815e..0708e769 100644
--- a/src/FSharp/Matrix.fs
+++ b/src/FSharp/Matrix.fs
@@ -135,14 +135,30 @@ module Matrix =
             for j=0 to M-1 do
                 C.[i,j] <- f i j (C.Item(i,j))
         C
-    
-    /// Map every matrix column into an element in a sequence using the given position dependent function.
-    let inline mapCols (f: int -> Vector<float> -> 'a) (A: #Matrix<float>) =
-        A.ColumnEnumerator() |> Seq.map (fun (j, col) -> f j col)
+        
+    /// In-place map every matrix column using the given position dependent function.
+    let inline inplaceMapCols (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 inplaceMapRows (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 column using the given position dependent function.
+    let inline mapCols (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        let A = A.Clone()
+        inplaceMapCols f A
+        A
 
-    /// Map every matrix row into an element in a sequence using the given position dependent function.
-    let inline mapRows (f: int -> Vector<float> -> 'a) (A: #Matrix<float>) =
-        A.RowEnumerator() |> Seq.map (fun (i, row) -> f i row)
+    /// 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()
+        inplaceMapRows f A
+        A
 
     /// In-place assignment.
     let inline inplaceAssign (f: int -> int -> float) (A: #Matrix<float>) =
@@ -170,13 +186,13 @@ module Matrix =
                 f <- f + A.Item(i,j)
         f
 
-    /// Returns the sum of the results generated by applying the function to each column of the matrix.
+    /// 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>) =
-        mapCols f A |> Seq.reduce (+)
+        A.ColumnEnumerator() |> Seq.map (fun (j,col) -> f j col) |> Seq.reduce (+)
     
-    /// Returns the sum of the results generated by applying the function to each row of the matrix.
+    /// 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>) =
-        mapRows f A |> Seq.reduce (+)
+        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>) =
diff --git a/src/FSharpUnitTests/Program.fs b/src/FSharpUnitTests/Program.fs
index 842be081..e15e6049 100644
--- a/src/FSharpUnitTests/Program.fs
+++ b/src/FSharpUnitTests/Program.fs
@@ -132,9 +132,9 @@ let MatrixTests =
         spec "Matrix.mapi"
             (Matrix.mapi (fun i j x -> float i * 100.0 + float j + x) largeM |> should equal (2.0 * largeM))
         spec "Matrix.mapCols"
-            (Matrix.mapCols (fun j col -> col.[0] + col.[1]) largeM |> Seq.toArray |> should array_equal [| for j in 0.0 .. 99.0 -> 2.0 * j + 100.0 |])
+            (Matrix.mapCols (fun j col -> col.Add(float j)) smallM |> should approximately_matrix_equal 14 (matrix [[0.3;1.3];[0.3;1.3]]))
         spec "Matrix.mapRows"
-            (Matrix.mapRows (fun i row -> row.[0] + row.[1]) largeM |> Seq.toArray |> should array_equal [| for i in 0.0 .. 99.0 -> 2.0 * 100.0 * i + 1.0 |])
+            (Matrix.mapRows (fun i row -> row.Add(float i)) smallM |> should approximately_matrix_equal 14 (matrix [[0.3;0.3];[1.3;1.3]]))
         spec "Matrix.inplaceAssign"
             ( let N = smallM.Clone()
               Matrix.inplaceAssign (fun i j -> 0.0) N