diff --git a/spatial/r3/mat.go b/spatial/r3/mat.go
index e638a1d8..eaa3f47f 100644
--- a/spatial/r3/mat.go
+++ b/spatial/r3/mat.go
@@ -198,6 +198,23 @@ func (m *Mat) Det() float64 {
 	return a*deta - b*detb + c*detc
 }
 
+// Skew sets the receiver to the 3×3 skew symmetric matrix
+// (right hand system) of v.
+//                  ⎡ 0 -z  y⎤
+//  Skew({x,y,z}) = ⎢ z  0 -x⎥
+//                  ⎣-y  x  0⎦
+func (m *Mat) Skew(v Vec) {
+	m.Set(0, 0, 0)
+	m.Set(0, 1, -v.Z)
+	m.Set(0, 2, v.Y)
+	m.Set(1, 0, v.Z)
+	m.Set(1, 1, 0)
+	m.Set(1, 2, -v.X)
+	m.Set(2, 0, -v.Y)
+	m.Set(2, 1, v.X)
+	m.Set(2, 2, 0)
+}
+
 // Hessian sets the receiver to the Hessian matrix of the scalar field at the point p,
 // approximated using finite differences with the given step sizes.
 // The field is evaluated at points in the area surrounding p by adding
diff --git a/spatial/r3/mat_safe.go b/spatial/r3/mat_safe.go
index 740259c1..1d23d2fe 100644
--- a/spatial/r3/mat_safe.go
+++ b/spatial/r3/mat_safe.go
@@ -60,6 +60,8 @@ func Eye() *Mat {
 //                  ⎡ 0 -z  y⎤
 //  Skew({x,y,z}) = ⎢ z  0 -x⎥
 //                  ⎣-y  x  0⎦
+//
+// DEPRECATED: use Mat.Skew()
 func Skew(v Vec) (M *Mat) {
 	return &Mat{&array{
 		0, -v.Z, v.Y,
diff --git a/spatial/r3/mat_test.go b/spatial/r3/mat_test.go
index 1ac3dd5f..d50d9c4d 100644
--- a/spatial/r3/mat_test.go
+++ b/spatial/r3/mat_test.go
@@ -103,9 +103,10 @@ func TestSkew(t *testing.T) {
 	const tol = 1e-16
 	rnd := rand.New(rand.NewSource(1))
 	for tc := 0; tc < 20; tc++ {
+		sk := NewMat(nil)
 		v1 := randomVec(rnd)
 		v2 := randomVec(rnd)
-		sk := Skew(v1)
+		sk.Skew(v1)
 		want := Cross(v1, v2)
 		got := sk.MulVec(v2)
 		if d := want.Sub(got); d.Dot(d) > tol {
diff --git a/spatial/r3/mat_unsafe.go b/spatial/r3/mat_unsafe.go
index c6c6d4f3..59b251f1 100644
--- a/spatial/r3/mat_unsafe.go
+++ b/spatial/r3/mat_unsafe.go
@@ -47,6 +47,8 @@ func Eye() *Mat {
 //                  ⎡ 0 -z  y⎤
 //  Skew({x,y,z}) = ⎢ z  0 -x⎥
 //                  ⎣-y  x  0⎦
+//
+// DEPRECATED: use Mat.Skew()
 func Skew(v Vec) (M *Mat) {
 	return &Mat{&array{
 		{0, -v.Z, v.Y},