spatial/r3: provide safe option for Mat operations

spatial/r3/mat.go

@@ -4,21 +4,12 @@
 
 package r3
 
-import (
+import "gonum.org/v1/gonum/mat"
-	"unsafe"
-
-	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/mat"
-)
-
-const (
-	badDim = "bad matrix dimensions"
-	badIdx = "bad matrix index"
-)
 
 // Mat represents a 3×3 matrix. Useful for rotation matrices and such.
+// The zero value is usable as the 3×3 zero matrix.
 type Mat struct {
-	data *[3][3]float64
+	data *array
 }
 
 var _ mat.Matrix = (*Mat)(nil)
@@ -27,58 +18,32 @@ var _ mat.Matrix = (*Mat)(nil)
 // with values passed as argument in row-major form. If val argument
 // is nil then NewMat returns a matrix filled with zeros.
 func NewMat(val []float64) *Mat {
-	if len(val) != 9 {
+	if len(val) == 9 {
+		return &Mat{arrayFrom(val)}
+	}
 	if val == nil {
-			return &Mat{data: new([3][3]float64)}
+		return &Mat{new(array)}
 	}
-		panic(badDim)
+	panic(mat.ErrShape)
-	}
-	m := Mat{}
-	m.setBackingSlice(val)
-	return &m
 }
 
 // Dims returns the number of rows and columns of this matrix.
 // This method will always return 3×3 for a Mat.
 func (m *Mat) Dims() (r, c int) { return 3, 3 }
 
-// At returns the value of a matrix element at row i, column j.
-// At expects indices in the range [0,2].
-// It will panic if i or j are out of bounds for the matrix.
-func (m *Mat) At(i, j int) float64 {
-	return m.data[i][j]
-}
-
-// Set sets the element at row i, column j to the value v.
-func (m *Mat) Set(i, j int, v float64) {
-	m.data[i][j] = v
-}
-
 // T returns the transpose of Mat. Changes in the receiver will be reflected in the returned matrix.
 func (m *Mat) T() mat.Matrix { return mat.Transpose{Matrix: m} }
 
-// RawMatrix returns the blas representation of the matrix with the backing data of this matrix.
-// Changes to the returned matrix will be reflected in the receiver.
-func (m *Mat) RawMatrix() blas64.General {
-	return blas64.General{Rows: 3, Cols: 3, Data: m.backingSlice(), Stride: 3}
-}
-
-// Eye returns the 3×3 Identity matrix
-func Eye() *Mat {
-	return &Mat{data: &[3][3]float64{
-		{1, 0, 0},
-		{0, 1, 0},
-		{0, 0, 1},
-	}}
-}
-
 // Scale multiplies the elements of a by f, placing the result in the receiver.
 //
 // See the mat.Scaler interface for more information.
 func (m *Mat) Scale(f float64, a mat.Matrix) {
 	r, c := a.Dims()
 	if r != 3 || c != 3 {
-		panic(badDim)
+		panic(mat.ErrShape)
+	}
+	if m.data == nil {
+		m.data = new(array)
 	}
 	for i := 0; i < 3; i++ {
 		for j := 0; j < 3; j++ {
@@ -90,6 +55,9 @@ func (m *Mat) Scale(f float64, a mat.Matrix) {
 // Performs matrix multiplication on v:
 //  result = M * v
 func (m *Mat) MulVec(v Vec) Vec {
+	if m.data == nil {
+		return Vec{}
+	}
 	return Vec{
 		X: v.X*m.At(0, 0) + v.Y*m.At(0, 1) + v.Z*m.At(0, 2),
 		Y: v.X*m.At(1, 0) + v.Y*m.At(1, 1) + v.Z*m.At(1, 2),
@@ -100,6 +68,9 @@ func (m *Mat) MulVec(v Vec) Vec {
 // Performs transposed matrix multiplication on v:
 //  result = Mᵀ * v
 func (m *Mat) MulVecTrans(v Vec) Vec {
+	if m.data == nil {
+		return Vec{}
+	}
 	return Vec{
 		X: v.X*m.At(0, 0) + v.Y*m.At(1, 0) + v.Z*m.At(2, 0),
 		Y: v.X*m.At(0, 1) + v.Y*m.At(1, 1) + v.Z*m.At(2, 1),
@@ -107,73 +78,15 @@ func (m *Mat) MulVecTrans(v Vec) Vec {
 	}
 }
 
-// Skew returns 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 Skew(v Vec) (M *Mat) {
-	return &Mat{data: &[3][3]float64{
-		{0, -v.Z, v.Y},
-		{v.Z, 0, -v.X},
-		{-v.Y, v.X, 0},
-	}}
-}
-
-// Mul takes the matrix product of a and b, placing the result in the receiver.
-// If the number of columns in a does not equal 3, Mul will panic.
-func (m *Mat) Mul(a, b mat.Matrix) {
-	ra, ca := a.Dims()
-	rb, cb := b.Dims()
-	switch {
-	case ra != 3:
-		panic(badDim)
-	case cb != 3:
-		panic(badDim)
-	case ca != rb:
-		panic(badDim)
-	}
-	if ca != 3 {
-		// General matrix multiplication for the case where the inner dimension is not 3.
-		t := mat.NewDense(3, 3, m.backingSlice())
-		t.Mul(a, b)
-		return
-	}
-
-	a00 := a.At(0, 0)
-	b00 := b.At(0, 0)
-	a01 := a.At(0, 1)
-	b01 := b.At(0, 1)
-	a02 := a.At(0, 2)
-	b02 := b.At(0, 2)
-	a10 := a.At(1, 0)
-	b10 := b.At(1, 0)
-	a11 := a.At(1, 1)
-	b11 := b.At(1, 1)
-	a12 := a.At(1, 2)
-	b12 := b.At(1, 2)
-	a20 := a.At(2, 0)
-	b20 := b.At(2, 0)
-	a21 := a.At(2, 1)
-	b21 := b.At(2, 1)
-	a22 := a.At(2, 2)
-	b22 := b.At(2, 2)
-	m.data[0][0] = a00*b00 + a01*b10 + a02*b20
-	m.data[0][1] = a00*b01 + a01*b11 + a02*b21
-	m.data[0][2] = a00*b02 + a01*b12 + a02*b22
-	m.data[1][0] = a10*b00 + a11*b10 + a12*b20
-	m.data[1][1] = a10*b01 + a11*b11 + a12*b21
-	m.data[1][2] = a10*b02 + a11*b12 + a12*b22
-	m.data[2][0] = a20*b00 + a21*b10 + a22*b20
-	m.data[2][1] = a20*b01 + a21*b11 + a22*b21
-	m.data[2][2] = a20*b02 + a21*b12 + a22*b22
-}
-
 // CloneFrom makes a copy of a into the receiver m.
 // Mat expects a 3×3 input matrix.
 func (m *Mat) CloneFrom(a mat.Matrix) {
 	r, c := a.Dims()
 	if r != 3 || c != 3 {
-		panic(badDim)
+		panic(mat.ErrShape)
+	}
+	if m.data == nil {
+		m.data = new(array)
 	}
 	for i := 0; i < 3; i++ {
 		for j := 0; j < 3; j++ {
@@ -186,47 +99,56 @@ func (m *Mat) CloneFrom(a mat.Matrix) {
 // Sub will panic if the two matrices do not have the same shape.
 func (m *Mat) Sub(a, b mat.Matrix) {
 	if r, c := a.Dims(); r != 3 || c != 3 {
-		panic(badDim)
+		panic(mat.ErrShape)
 	}
 	if r, c := b.Dims(); r != 3 || c != 3 {
-		panic(badDim)
+		panic(mat.ErrShape)
+	}
+	if m.data == nil {
+		m.data = new(array)
 	}
 
-	m.data[0][0] = a.At(0, 0) - b.At(0, 0)
+	m.Set(0, 0, a.At(0, 0)-b.At(0, 0))
-	m.data[0][1] = a.At(0, 1) - b.At(0, 1)
+	m.Set(0, 1, a.At(0, 1)-b.At(0, 1))
-	m.data[0][2] = a.At(0, 2) - b.At(0, 2)
+	m.Set(0, 2, a.At(0, 2)-b.At(0, 2))
-	m.data[1][0] = a.At(1, 0) - b.At(1, 0)
+	m.Set(1, 0, a.At(1, 0)-b.At(1, 0))
-	m.data[1][1] = a.At(1, 1) - b.At(1, 1)
+	m.Set(1, 1, a.At(1, 1)-b.At(1, 1))
-	m.data[1][2] = a.At(1, 2) - b.At(1, 2)
+	m.Set(1, 2, a.At(1, 2)-b.At(1, 2))
-	m.data[2][0] = a.At(2, 0) - b.At(2, 0)
+	m.Set(2, 0, a.At(2, 0)-b.At(2, 0))
-	m.data[2][1] = a.At(2, 1) - b.At(2, 1)
+	m.Set(2, 1, a.At(2, 1)-b.At(2, 1))
-	m.data[2][2] = a.At(2, 2) - b.At(2, 2)
+	m.Set(2, 2, a.At(2, 2)-b.At(2, 2))
 }
 
 // Add adds a and b element-wise, placing the result in the receiver. Add will panic if the two matrices do not have the same shape.
 func (m *Mat) Add(a, b mat.Matrix) {
 	if r, c := a.Dims(); r != 3 || c != 3 {
-		panic(badDim)
+		panic(mat.ErrShape)
 	}
 	if r, c := b.Dims(); r != 3 || c != 3 {
-		panic(badDim)
+		panic(mat.ErrShape)
+	}
+	if m.data == nil {
+		m.data = new(array)
 	}
 
-	m.data[0][0] = a.At(0, 0) + b.At(0, 0)
+	m.Set(0, 0, a.At(0, 0)+b.At(0, 0))
-	m.data[0][1] = a.At(0, 1) + b.At(0, 1)
+	m.Set(0, 1, a.At(0, 1)+b.At(0, 1))
-	m.data[0][2] = a.At(0, 2) + b.At(0, 2)
+	m.Set(0, 2, a.At(0, 2)+b.At(0, 2))
-	m.data[1][0] = a.At(1, 0) + b.At(1, 0)
+	m.Set(1, 0, a.At(1, 0)+b.At(1, 0))
-	m.data[1][1] = a.At(1, 1) + b.At(1, 1)
+	m.Set(1, 1, a.At(1, 1)+b.At(1, 1))
-	m.data[1][2] = a.At(1, 2) + b.At(1, 2)
+	m.Set(1, 2, a.At(1, 2)+b.At(1, 2))
-	m.data[2][0] = a.At(2, 0) + b.At(2, 0)
+	m.Set(2, 0, a.At(2, 0)+b.At(2, 0))
-	m.data[2][1] = a.At(2, 1) + b.At(2, 1)
+	m.Set(2, 1, a.At(2, 1)+b.At(2, 1))
-	m.data[2][2] = a.At(2, 2) + b.At(2, 2)
+	m.Set(2, 2, a.At(2, 2)+b.At(2, 2))
 }
 
 // VecRow returns the elements in the ith row of the receiver.
 func (m *Mat) VecRow(i int) Vec {
 	if i > 2 {
-		panic(badIdx)
+		panic(mat.ErrRowAccess)
+	}
+	if m.data == nil {
+		return Vec{}
 	}
 	return Vec{X: m.At(i, 0), Y: m.At(i, 1), Z: m.At(i, 2)}
 }
@@ -234,17 +156,10 @@ func (m *Mat) VecRow(i int) Vec {
 // VecCol returns the elements in the jth column of the receiver.
 func (m *Mat) VecCol(j int) Vec {
 	if j > 2 {
-		panic(badIdx)
+		panic(mat.ErrColAccess)
+	}
+	if m.data == nil {
+		return Vec{}
 	}
 	return Vec{X: m.At(0, j), Y: m.At(1, j), Z: m.At(2, j)}
 }
-
-// setBackingSlice requires unsafe.
-func (m *Mat) setBackingSlice(vals []float64) {
-	m.data = (*[3][3]float64)(unsafe.Pointer(&vals[0]))
-}
-
-// backingSlice requires unsafe.
-func (m *Mat) backingSlice() []float64 {
-	return (*[9]float64)(unsafe.Pointer(m.data))[:]
-}

spatial/r3/mat_safe.go

@@ -0,0 +1,138 @@
+// Copyright ©2021 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+//go:build safe
+// +build safe
+
+// TODO(kortschak): Get rid of this rigmarole if https://golang.org/issue/50118
+// is accepted.
+
+package r3
+
+import (
+	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/mat"
+)
+
+type array [9]float64
+
+// At returns the value of a matrix element at row i, column j.
+// At expects indices in the range [0,2].
+// It will panic if i or j are out of bounds for the matrix.
+func (m *Mat) At(i, j int) float64 {
+	if uint(i) > 2 {
+		panic(mat.ErrRowAccess)
+	}
+	if uint(j) > 2 {
+		panic(mat.ErrColAccess)
+	}
+	if m.data == nil {
+		m.data = new(array)
+	}
+	return m.data[i*3+j]
+}
+
+// Set sets the element at row i, column j to the value v.
+func (m *Mat) Set(i, j int, v float64) {
+	if uint(i) > 2 {
+		panic(mat.ErrRowAccess)
+	}
+	if uint(j) > 2 {
+		panic(mat.ErrColAccess)
+	}
+	if m.data == nil {
+		m.data = new(array)
+	}
+	m.data[i*3+j] = v
+}
+
+// Eye returns the 3×3 Identity matrix
+func Eye() *Mat {
+	return &Mat{&array{
+		1, 0, 0,
+		0, 1, 0,
+		0, 0, 1,
+	}}
+}
+
+// Skew returns 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 Skew(v Vec) (M *Mat) {
+	return &Mat{&array{
+		0, -v.Z, v.Y,
+		v.Z, 0, -v.X,
+		-v.Y, v.X, 0,
+	}}
+}
+
+// Mul takes the matrix product of a and b, placing the result in the receiver.
+// If the number of columns in a does not equal 3, Mul will panic.
+func (m *Mat) Mul(a, b mat.Matrix) {
+	ra, ca := a.Dims()
+	rb, cb := b.Dims()
+	switch {
+	case ra != 3:
+		panic(mat.ErrShape)
+	case cb != 3:
+		panic(mat.ErrShape)
+	case ca != rb:
+		panic(mat.ErrShape)
+	}
+	if m.data == nil {
+		m.data = new(array)
+	}
+	if ca != 3 {
+		// General matrix multiplication for the case where the inner dimension is not 3.
+		var t mat.Dense
+		t.Mul(a, b)
+		copy(m.data[:], t.RawMatrix().Data)
+		return
+	}
+
+	a00 := a.At(0, 0)
+	b00 := b.At(0, 0)
+	a01 := a.At(0, 1)
+	b01 := b.At(0, 1)
+	a02 := a.At(0, 2)
+	b02 := b.At(0, 2)
+	a10 := a.At(1, 0)
+	b10 := b.At(1, 0)
+	a11 := a.At(1, 1)
+	b11 := b.At(1, 1)
+	a12 := a.At(1, 2)
+	b12 := b.At(1, 2)
+	a20 := a.At(2, 0)
+	b20 := b.At(2, 0)
+	a21 := a.At(2, 1)
+	b21 := b.At(2, 1)
+	a22 := a.At(2, 2)
+	b22 := b.At(2, 2)
+	*(m.data) = array{
+		a00*b00 + a01*b10 + a02*b20, a00*b01 + a01*b11 + a02*b21, a00*b02 + a01*b12 + a02*b22,
+		a10*b00 + a11*b10 + a12*b20, a10*b01 + a11*b11 + a12*b21, a10*b02 + a11*b12 + a12*b22,
+		a20*b00 + a21*b10 + a22*b20, a20*b01 + a21*b11 + a22*b21, a20*b02 + a21*b12 + a22*b22,
+	}
+}
+
+// RawMatrix returns the blas representation of the matrix with the backing
+// data of this matrix. Changes to the returned matrix will be reflected in
+// the receiver.
+func (m *Mat) RawMatrix() blas64.General {
+	if m.data == nil {
+		m.data = new(array)
+	}
+	return blas64.General{Rows: 3, Cols: 3, Data: m.data[:], Stride: 3}
+}
+
+// BUG(kortschak): Implementing NewMat without unsafe conversion or slice to
+// array pointer conversion leaves it with semantics that do not match the
+// sharing semantics that exist in the mat package.
+func arrayFrom(vals []float64) *array {
+	// TODO(kortschak): Use array conversion when go1.16 is no longer supported.
+	var a array
+	copy(a[:], vals)
+	return &a
+}

spatial/r3/mat_test.go

@@ -5,7 +5,6 @@
 package r3
 
 import (
-	"math"
 	"testing"
 
 	"golang.org/x/exp/rand"
@@ -13,55 +12,108 @@ import (
 	"gonum.org/v1/gonum/mat"
 )
 
+func TestMatAdd(t *testing.T) {
+	const tol = 1e-16
+	rnd := rand.New(rand.NewSource(1))
+	for tc := 0; tc < 20; tc++ {
+		a := randomMat(rnd)
+		b := randomMat(rnd)
+		var (
+			want mat.Dense
+			got  Mat
+		)
+		want.Add(a, b)
+		got.Add(a, b)
+		if !mat.EqualApprox(&got, &want, tol) {
+			t.Errorf("unexpected result for matrix add:\ngot:\n%v\nwant:\n%v", mat.Formatted(&got), mat.Formatted(&want))
+		}
+	}
+}
+
+func TestMatSub(t *testing.T) {
+	const tol = 1e-16
+	rnd := rand.New(rand.NewSource(1))
+	for tc := 0; tc < 20; tc++ {
+		a := randomMat(rnd)
+		b := randomMat(rnd)
+		var (
+			want mat.Dense
+			got  Mat
+		)
+		want.Sub(a, b)
+		got.Sub(a, b)
+		if !mat.EqualApprox(&got, &want, tol) {
+			t.Errorf("unexpected result for matrix subtract:\ngot:\n%v\nwant:\n%v", mat.Formatted(&got), mat.Formatted(&want))
+		}
+	}
+}
+
+func TestMatMul(t *testing.T) {
+	const tol = 1e-14
+	rnd := rand.New(rand.NewSource(1))
+	for tc := 0; tc < 20; tc++ {
+		a := randomMat(rnd)
+		b := randomMat(rnd)
+		var (
+			want mat.Dense
+			got  Mat
+		)
+		want.Mul(a, b)
+		got.Mul(a, b)
+		if !mat.EqualApprox(&got, &want, tol) {
+			t.Errorf("unexpected result for matrix multiply:\ngot:\n%v\nwant:\n%v", mat.Formatted(&got), mat.Formatted(&want))
+		}
+	}
+}
+
 func TestMatScale(t *testing.T) {
-	const tol = 1e-12
+	const tol = 1e-16
 	rnd := rand.New(rand.NewSource(1))
 	for tc := 0; tc < 20; tc++ {
 		v := rnd.Float64()
 		a := randomMat(rnd)
-		gotmat := NewMat(nil)
+		var (
-		gotmat.Scale(v, a)
+			want mat.Dense
-		for iv := range a.data {
+			got  Mat
-			i := iv / 3
+		)
-			j := iv % 3
+		want.Scale(v, a)
-			expect := v * a.At(i, j)
+		got.Scale(v, a)
-			got := gotmat.At(i, j)
+		if !mat.EqualApprox(&got, &want, tol) {
-			if math.Abs(got-expect) > tol {
+			t.Errorf("unexpected result for matrix scale:\ngot:\n%v\nwant:\n%v", mat.Formatted(&got), mat.Formatted(&want))
-				t.Errorf(
-					"case %d: got=%v, want=%v",
-					tc, got, expect)
-			}
 		}
 	}
 }
 
 func TestMatCloneFrom(t *testing.T) {
+	const tol = 1e-16
 	rnd := rand.New(rand.NewSource(1))
 	for tc := 0; tc < 20; tc++ {
-		a := randomMat(rnd)
+		want := randomMat(rnd)
-		gotmat := NewMat(nil)
+		got := NewMat(nil)
-		gotmat.CloneFrom(a)
+		got.CloneFrom(want)
-		if !mat.Equal(a, gotmat) {
+		if !mat.EqualApprox(got, want, tol) {
-			t.Error("Clonefrom fail")
+			t.Errorf("unexpected result from CloneFrom:\ngot:\n%v\nwant:\n%v", mat.Formatted(got), mat.Formatted(want))
 		}
 	}
 }
 
 func TestSkew(t *testing.T) {
+	const tol = 1e-16
 	rnd := rand.New(rand.NewSource(1))
 	for tc := 0; tc < 20; tc++ {
 		v1 := randomVec(rnd)
 		v2 := randomVec(rnd)
 		sk := Skew(v1)
+		want := Cross(v1, v2)
 		got := sk.MulVec(v2)
-		expect := Cross(v1, v2)
+		if d := want.Sub(got); d.Dot(d) > tol {
-		if got != expect {
+			t.Errorf("r3.Cross(v1,v2) does not agree with r3.Skew(v1)*v2: got:%v want:%v", got, want)
-			t.Error("r3.Cross(v1,v2) not match with r3.Skew(v1)*v2")
 		}
 	}
 }
 
 func TestTranspose(t *testing.T) {
+	const tol = 1e-16
 	rnd := rand.New(rand.NewSource(1))
 	for tc := 0; tc < 20; tc++ {
 		d := mat.NewDense(3, 3, nil)
@@ -70,23 +122,24 @@ func TestTranspose(t *testing.T) {
 		mt := m.T()
 		dt := d.T()
 		if !mat.Equal(mt, dt) {
-			t.Error("Dense.T() not equal to r3.Mat.T()")
+			t.Errorf("Dense.T() not equal to r3.Mat.T():\ngot:\n%v\nwant:\n%v", mat.Formatted(mt), mat.Formatted(dt))
 		}
 		vd := mat.NewVecDense(3, nil)
 		v := randomVec(rnd)
 		vd.SetVec(0, v.X)
 		vd.SetVec(1, v.Y)
 		vd.SetVec(2, v.Z)
-		got := m.MulVecTrans(v)
 		vd.MulVec(dt, vd)
-		if vd.AtVec(0) != got.X || vd.AtVec(1) != got.Y || vd.AtVec(2) != got.Z {
+		want := Vec{X: vd.AtVec(0), Y: vd.AtVec(1), Z: vd.AtVec(2)}
-			t.Error("VecDense.MulVec(dense.T()) not equal to r3.Mat.MulVec(r3.Vec)")
+		got := m.MulVecTrans(v)
+		if d := want.Sub(got); d.Dot(d) > tol {
+			t.Errorf("VecDense.MulVec(dense.T()) not agree with r3.Mat.MulVec(r3.Vec): got:%v want:%v", got, want)
 		}
 	}
 }
 
 func randomMat(rnd *rand.Rand) *Mat {
-	m := Mat{data: new([3][3]float64)}
+	m := Mat{new(array)}
 	for iv := 0; iv < 9; iv++ {
 		i := iv / 3
 		j := iv % 3

spatial/r3/mat_unsafe.go

@@ -0,0 +1,127 @@
+// Copyright ©2021 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+//go:build !safe
+// +build !safe
+
+package r3
+
+import (
+	"unsafe"
+
+	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/mat"
+)
+
+type array [3][3]float64
+
+// At returns the value of a matrix element at row i, column j.
+// At expects indices in the range [0,2].
+// It will panic if i or j are out of bounds for the matrix.
+func (m *Mat) At(i, j int) float64 {
+	if m.data == nil {
+		m.data = new(array)
+	}
+	return m.data[i][j]
+}
+
+// Set sets the element at row i, column j to the value v.
+func (m *Mat) Set(i, j int, v float64) {
+	if m.data == nil {
+		m.data = new(array)
+	}
+	m.data[i][j] = v
+}
+
+// Eye returns the 3×3 Identity matrix
+func Eye() *Mat {
+	return &Mat{&array{
+		{1, 0, 0},
+		{0, 1, 0},
+		{0, 0, 1},
+	}}
+}
+
+// Skew returns 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 Skew(v Vec) (M *Mat) {
+	return &Mat{&array{
+		{0, -v.Z, v.Y},
+		{v.Z, 0, -v.X},
+		{-v.Y, v.X, 0},
+	}}
+}
+
+// Mul takes the matrix product of a and b, placing the result in the receiver.
+// If the number of columns in a does not equal 3, Mul will panic.
+func (m *Mat) Mul(a, b mat.Matrix) {
+	ra, ca := a.Dims()
+	rb, cb := b.Dims()
+	switch {
+	case ra != 3:
+		panic(mat.ErrShape)
+	case cb != 3:
+		panic(mat.ErrShape)
+	case ca != rb:
+		panic(mat.ErrShape)
+	}
+	if m.data == nil {
+		m.data = new(array)
+	}
+	if ca != 3 {
+		// General matrix multiplication for the case where the inner dimension is not 3.
+		t := mat.NewDense(3, 3, m.slice())
+		t.Mul(a, b)
+		return
+	}
+
+	a00 := a.At(0, 0)
+	b00 := b.At(0, 0)
+	a01 := a.At(0, 1)
+	b01 := b.At(0, 1)
+	a02 := a.At(0, 2)
+	b02 := b.At(0, 2)
+	a10 := a.At(1, 0)
+	b10 := b.At(1, 0)
+	a11 := a.At(1, 1)
+	b11 := b.At(1, 1)
+	a12 := a.At(1, 2)
+	b12 := b.At(1, 2)
+	a20 := a.At(2, 0)
+	b20 := b.At(2, 0)
+	a21 := a.At(2, 1)
+	b21 := b.At(2, 1)
+	a22 := a.At(2, 2)
+	b22 := b.At(2, 2)
+	m.data[0][0] = a00*b00 + a01*b10 + a02*b20
+	m.data[0][1] = a00*b01 + a01*b11 + a02*b21
+	m.data[0][2] = a00*b02 + a01*b12 + a02*b22
+	m.data[1][0] = a10*b00 + a11*b10 + a12*b20
+	m.data[1][1] = a10*b01 + a11*b11 + a12*b21
+	m.data[1][2] = a10*b02 + a11*b12 + a12*b22
+	m.data[2][0] = a20*b00 + a21*b10 + a22*b20
+	m.data[2][1] = a20*b01 + a21*b11 + a22*b21
+	m.data[2][2] = a20*b02 + a21*b12 + a22*b22
+}
+
+// RawMatrix returns the blas representation of the matrix with the backing
+// data of this matrix. Changes to the returned matrix will be reflected in
+// the receiver.
+func (m *Mat) RawMatrix() blas64.General {
+	if m.data == nil {
+		m.data = new(array)
+	}
+	return blas64.General{Rows: 3, Cols: 3, Data: m.slice(), Stride: 3}
+}
+
+func arrayFrom(vals []float64) *array {
+	// TODO(kortschak): Use array conversion when go1.16 is no longer supported.
+	return (*array)(unsafe.Pointer(&vals[0]))
+}
+
+func (m *Mat) slice() []float64 {
+	return (*[9]float64)(unsafe.Pointer(m.data))[:]
+}