Repository: gonum/gonum (language: go)

Issue:
spatial/r3: Adding a `r3.Mat` type
I work with rotation tensors and vectors of dimension 3. I found myself writing a package for matrix and r3.Vec interoperability.



<!--
Please make sure your issue title matches the Go convention; a summary
of the problem, prefixed by "proposal:" and the primary affected package.

It may be worth discussing the proposal at the Gonum discussion list
before filing the proposal issue:

https://groups.google.com/forum/#!forum/gonum-dev
-->
### Background
I work with rotation tensors and vectors of dimension 3. I found myself writing a package for matrix and r3.Vec interoperability.
These tensors are used for attitude representation of a launch vehicle and are the industry standard way of representing this concept. One can use euler angles or quaternions to construct/decompose these matrices. See _Modeling and simulation of aerospace vehicle dynamics_ by P. Zipfel.

### Proposal
I propose the following addition to `spatial/r3`

<details>

```go
import (
	"fmt"
	"math"

	"gonum.org/v1/gonum/mat"
	"gonum.org/v1/gonum/spatial/r3"
)

const (
	badDim = "bad matrix dimensions"
	badIdx = "bad matrix index"
)

type Mat3 [3 * 3]float64

func NewMat3(val []float64) *Mat3 {
	if len(val) != 9 {
		if val == nil {
			return &Mat3{}
		}
		panic(badDim)
	}
	M := Mat3{}
	copy(M[:], val)
	return &M
}

// IdentityMat3 returns the Identity 3x3 matrix
func IdentityMat3() *Mat3 { return &Mat3{1, 0, 0, 0, 1, 0, 0, 0, 1} }

func (m Mat3) String() string {
	for i := range m {
		if math.Abs(m[i]) < 1e-6 {
			m[i] = 0
		}
	}
	return fmt.Sprintf("[%0.6g\t%0.6g\t%0.6g]\n[%0.6g\t%0.6g\t%0.6g]\n[%0.6g\t%0.6g\t%0.6g]",
		m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8])
}

func (m Mat3) CloneDense() *mat.Dense {
	M := mat.NewDense(3, 3, m[:])
	return M
	// for i := 0; i < 3; i++ {
	// 	for j := 0; j < 3; j++ {
	// 		M.Set(i, j, m.at(i, j))
	// 	}
	// }
	// return M
}

// CopyToMat3 copies elements from src Matrix to
// dst Mat3.
func CopyToMat3(dst *Mat3, src mat.Matrix) {
	r, c := src.Dims()
	if r != 3 || c != 3 {
		panic(badDim)
	}
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			dst.set(i, j, src.At(i, j))
		}
	}
}

func (m *Mat3) Dims() (r, c int) { return 3, 3 }

func (m *Mat3) At(i, j int) float64 {
	if uint(i) > 2 || uint(j) > 2 {
		panic(badIdx)
	}
	return m.at(i, j)
}

func (m *Mat3) at(i, j int) float64 {
	return m[i*3+j]
}

func (m *Mat3) T() mat.Matrix { return mat.Transpose{Matrix: m} }
func (m *Mat3) Set(i, j int, v float64) {
	if uint(i) > 2 || uint(j) > 2 {
		panic(badIdx)
	}
	m.set(i, j, v)
}

// set sets element at position i,j without bounds checking
func (m *Mat3) set(i, j int, v float64) {
	m[i*3+j] = v
}

// Scale multiplies the elements of a by f, placing the result in the receiver.
//
// See the Scaler interface for more information.
func (m *Mat3) Scale(f float64, a mat.Matrix) {
	r, c := a.Dims()
	if r != 3 || c != 3 {
		panic(badDim)
	}
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			m.set(i, j, f*a.At(i, j))
		}
	}
}

// Performs matrix multiplication on v:
//  result = M * v
func (m *Mat3) MulVec(v r3.Vec) (result r3.Vec) {
	result.X = v.X*m[0] + v.Y*m[1] + v.Z*m[2]
	result.Y = v.X*m[3] + v.Y*m[4] + v.Z*m[5]
	result.Z = v.X*m[6] + v.Y*m[7] + v.Z*m[8]
	return result
}

// Performs transposed matrix multiplication on v:
//  result = Mᵀ * v
func (m *Mat3) MulVecTrans(v r3.Vec) (result r3.Vec) {
	result.X = v.X*m[0] + v.Y*m[3] + v.Z*m[6]
	result.Y = v.X*m[1] + v.Y*m[4] + v.Z*m[7]
	result.Z = v.X*m[2] + v.Y*m[5] + v.Z*m[8]
	return result
}

// Skew returns the skew symmetric matrix (right hand system) of v.
//                  [0  -z  y]
//  Skew({x,y,z}) = [z   0 -x]
//                  [-y  x  0]
func Skew(v r3.Vec) (M *Mat3) {
	return &Mat3{
		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 *Mat3) Mul(a, b mat.Matrix) {
	// TODO: add optimized operations for Mat3*Mat3
	// A, oka := a.(*Mat3)
	// B, okb := b.(*Mat3)
	// if oka && okb {
	// 	m.mul(A, B)
	// 	return
	// }
	mat.NewDense(3, 3, m[:]).Mul(a, b)
}

// Subtracts b from a and stores result in receiver. Sub will panic if the two matrices do not have the same shape.
func (m *Mat3) Sub(a, b mat.Matrix) {
	mat.NewDense(3, 3, m[:]).Sub(a, b)
}

// 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 *Mat3) Add(a, b mat.Matrix) {
	mat.NewDense(3, 3, m[:]).Add(a, b)
}

func (m *Mat3) VecRow(i int) r3.Vec {
	if uint(i) > 2 {
		panic(badIdx)
	}
	return r3.Vec{X: m.at(i, 0), Y: m.at(i, 1), Z: m.at(i, 2)}
}

func (m *Mat3) VecCol(j int) r3.Vec {
	if uint(j) > 2 {
		panic(badIdx)
	}
	return r3.Vec{X: m.at(0, j), Y: m.at(1, j), Z: m.at(2, j)}
}
```

</details>

### Potential impact of proposal
<!-- Extent of change, amount of work to implement, who will do the work? etc. -->

Above is a working example of a 3x3 matrix. Of course with time the routines could be optimized for the special case of 3x3 arguments to Add, Mul, etc. Let me know what you think!

`Mat` implements the `mat.Matrix` interface. Maybe it is worth having `r3.Vec` implement the `Vector` interface so that there is better gonum package interoperability.


Discussion / hints:
I am in favor of this proposal (I vaguely remember making the point for such a thing (well, `Vec` actually) last year).

a couple of "drive-by" comments:
- it would be great to be able to implement `mat.Matrix` w/o importing `mat` (but that's not possible)
- I would implement `Mat3` as a `struct{ x00, x01, ... x22 float64}` as it's more amenable to SSA optimizations
- I would `s/Mat3/Mat3x3/g`
- I would `s/IdentityMat3/Id3x3/` or `Identity3x3` or `Eye3x3`
- I would drop `CloneDense` (that can be done with `mat.DenseCopyFrom`)
- I would drop `CopyToMat3` and instead implement the `mat.ClonerFrom` interface

@kortschak WDYT?
I think we can probably simplify the name to `Mat` since r3 is already in the package name. With that in mind, `IdentityMat` and I agree with the drops.

I am not convinced by the use of a struct with named fields for elements, I was thinking `struct { data [3][3]float64 }` all on pointers. There is no need for `String`
I think changing `type Mat [3*3]float64` for `type Mat struct {data [3][3]float64}` doesn't buy you anything, performance-wise.
at least that's what I understood from this discussion:
https://github.com/gonum/gonum/pull/1297#issuecomment-603466614
Yeah, I wasn't thinking about performance. It is clearer though. Thanks for linking that, I thought I remembered a conversation about this kind of thing. If it is hidden behind a struct wall, it can be changed.
@sbinet @kortschak  I agree with all your comments. I can make a PR for the agreed upon drops. 

Fix patch:
diff --git a/spatial/r3/mat.go b/spatial/r3/mat.go
new file mode 100644
index 000000000..c2e991155
--- /dev/null
+++ b/spatial/r3/mat.go
@@ -0,0 +1,250 @@
+// 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.
+
+package r3
+
+import (
+	"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.
+type Mat struct {
+	data *[3][3]float64
+}
+
+var _ mat.Matrix = (*Mat)(nil)
+
+// NewMat returns a new 3×3 matrix Mat type and populates its elements
+// 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 val == nil {
+			return &Mat{data: new([3][3]float64)}
+		}
+		panic(badDim)
+	}
+	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)
+	}
+	for i := 0; i < 3; i++ {
+		for j := 0; j < 3; j++ {
+			m.Set(i, j, f*a.At(i, j))
+		}
+	}
+}
+
+// Performs matrix multiplication on v:
+//  result = M * v
+func (m *Mat) MulVec(v Vec) 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),
+		Z: v.X*m.At(2, 0) + v.Y*m.At(2, 1) + v.Z*m.At(2, 2),
+	}
+}
+
+// Performs transposed matrix multiplication on v:
+//  result = Mᵀ * v
+func (m *Mat) MulVecTrans(v Vec) 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),
+		Z: v.X*m.At(0, 2) + v.Y*m.At(1, 2) + v.Z*m.At(2, 2),
+	}
+}
+
+// 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)
+	}
+	for i := 0; i < 3; i++ {
+		for j := 0; j < 3; j++ {
+			m.Set(i, j, a.At(i, j))
+		}
+	}
+}
+
+// Sub subtracts the matrix b from a, placing the result in the receiver.
+// 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)
+	}
+	if r, c := b.Dims(); r != 3 || c != 3 {
+		panic(badDim)
+	}
+
+	m.data[0][0] = a.At(0, 0) - b.At(0, 0)
+	m.data[0][1] = a.At(0, 1) - b.At(0, 1)
+	m.data[0][2] = a.At(0, 2) - b.At(0, 2)
+	m.data[1][0] = a.At(1, 0) - b.At(1, 0)
+	m.data[1][1] = a.At(1, 1) - b.At(1, 1)
+	m.data[1][2] = a.At(1, 2) - b.At(1, 2)
+	m.data[2][0] = a.At(2, 0) - b.At(2, 0)
+	m.data[2][1] = a.At(2, 1) - b.At(2, 1)
+	m.data[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)
+	}
+	if r, c := b.Dims(); r != 3 || c != 3 {
+		panic(badDim)
+	}
+
+	m.data[0][0] = a.At(0, 0) + b.At(0, 0)
+	m.data[0][1] = a.At(0, 1) + b.At(0, 1)
+	m.data[0][2] = a.At(0, 2) + b.At(0, 2)
+	m.data[1][0] = a.At(1, 0) + b.At(1, 0)
+	m.data[1][1] = a.At(1, 1) + b.At(1, 1)
+	m.data[1][2] = a.At(1, 2) + b.At(1, 2)
+	m.data[2][0] = a.At(2, 0) + b.At(2, 0)
+	m.data[2][1] = a.At(2, 1) + b.At(2, 1)
+	m.data[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)
+	}
+	return Vec{X: m.At(i, 0), Y: m.At(i, 1), Z: m.At(i, 2)}
+}
+
+// VecCol returns the elements in the jth column of the receiver.
+func (m *Mat) VecCol(j int) Vec {
+	if j > 2 {
+		panic(badIdx)
+	}
+	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))[:]
+}


Test patch:
diff --git a/spatial/r3/mat_test.go b/spatial/r3/mat_test.go
new file mode 100644
index 000000000..66e81330e
--- /dev/null
+++ b/spatial/r3/mat_test.go
@@ -0,0 +1,103 @@
+// 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.
+
+package r3
+
+import (
+	"math"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/mat"
+)
+
+func TestMatScale(t *testing.T) {
+	const tol = 1e-12
+	rnd := rand.New(rand.NewSource(1))
+	for tc := 0; tc < 20; tc++ {
+		v := rnd.Float64()
+		a := randomMat(rnd)
+		gotmat := NewMat(nil)
+		gotmat.Scale(v, a)
+		for iv := range a.data {
+			i := iv / 3
+			j := iv % 3
+			expect := v * a.At(i, j)
+			got := gotmat.At(i, j)
+			if math.Abs(got-expect) > tol {
+				t.Errorf(
+					"case %d: got=%v, want=%v",
+					tc, got, expect)
+			}
+		}
+	}
+}
+
+func TestMatCloneFrom(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for tc := 0; tc < 20; tc++ {
+		a := randomMat(rnd)
+		gotmat := NewMat(nil)
+		gotmat.CloneFrom(a)
+		if !mat.Equal(a, gotmat) {
+			t.Error("Clonefrom fail")
+		}
+	}
+}
+
+func TestSkew(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for tc := 0; tc < 20; tc++ {
+		v1 := randomVec(rnd)
+		v2 := randomVec(rnd)
+		sk := Skew(v1)
+		got := sk.MulVec(v2)
+		expect := Cross(v1, v2)
+		if got != expect {
+			t.Error("r3.Cross(v1,v2) not match with r3.Skew(v1)*v2")
+		}
+	}
+}
+
+func TestTranspose(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for tc := 0; tc < 20; tc++ {
+		d := mat.NewDense(3, 3, nil)
+		m := randomMat(rnd)
+		d.CloneFrom(m)
+		mt := m.T()
+		dt := d.T()
+		if !mat.Equal(mt, dt) {
+			t.Error("Dense.T() not equal to r3.Mat.T()")
+		}
+		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 {
+			t.Error("VecDense.MulVec(dense.T()) not equal to r3.Mat.MulVec(r3.Vec)")
+		}
+	}
+}
+
+func randomMat(rnd *rand.Rand) *Mat {
+	m := Mat{data: new([3][3]float64)}
+	for iv := 0; iv < 9; iv++ {
+		i := iv / 3
+		j := iv % 3
+		m.Set(i, j, (rnd.Float64()-0.5)*20)
+	}
+	return &m
+}
+
+func randomVec(rnd *rand.Rand) (v Vec) {
+	v.X = (rnd.Float64() - 0.5) * 20
+	v.Y = (rnd.Float64() - 0.5) * 20
+	v.Z = (rnd.Float64() - 0.5) * 20
+	return v
+}