diff --git a/spatial/r3/mat.go b/spatial/r3/mat.go index c2e99115..37acc155 100644 --- a/spatial/r3/mat.go +++ b/spatial/r3/mat.go @@ -4,21 +4,12 @@ package r3 -import ( - "unsafe" - - "gonum.org/v1/gonum/blas/blas64" - "gonum.org/v1/gonum/mat" -) - -const ( - badDim = "bad matrix dimensions" - badIdx = "bad matrix index" -) +import "gonum.org/v1/gonum/mat" // 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 val == nil { - return &Mat{data: new([3][3]float64)} - } - panic(badDim) + if len(val) == 9 { + return &Mat{arrayFrom(val)} } - m := Mat{} - m.setBackingSlice(val) - return &m + if val == nil { + return &Mat{new(array)} + } + panic(mat.ErrShape) } // 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.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) + m.Set(0, 0, a.At(0, 0)-b.At(0, 0)) + m.Set(0, 1, a.At(0, 1)-b.At(0, 1)) + m.Set(0, 2, a.At(0, 2)-b.At(0, 2)) + m.Set(1, 0, a.At(1, 0)-b.At(1, 0)) + m.Set(1, 1, a.At(1, 1)-b.At(1, 1)) + m.Set(1, 2, a.At(1, 2)-b.At(1, 2)) + m.Set(2, 0, a.At(2, 0)-b.At(2, 0)) + m.Set(2, 1, a.At(2, 1)-b.At(2, 1)) + 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.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) + m.Set(0, 0, a.At(0, 0)+b.At(0, 0)) + m.Set(0, 1, a.At(0, 1)+b.At(0, 1)) + m.Set(0, 2, a.At(0, 2)+b.At(0, 2)) + m.Set(1, 0, a.At(1, 0)+b.At(1, 0)) + m.Set(1, 1, a.At(1, 1)+b.At(1, 1)) + m.Set(1, 2, a.At(1, 2)+b.At(1, 2)) + m.Set(2, 0, a.At(2, 0)+b.At(2, 0)) + m.Set(2, 1, a.At(2, 1)+b.At(2, 1)) + 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))[:] -} diff --git a/spatial/r3/mat_safe.go b/spatial/r3/mat_safe.go new file mode 100644 index 00000000..ca3e7c5a --- /dev/null +++ b/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 +} diff --git a/spatial/r3/mat_test.go b/spatial/r3/mat_test.go index 66e81330..eb165b01 100644 --- a/spatial/r3/mat_test.go +++ b/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) - 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) - } + var ( + want mat.Dense + got Mat + ) + want.Scale(v, a) + got.Scale(v, a) + if !mat.EqualApprox(&got, &want, tol) { + t.Errorf("unexpected result for matrix scale:\ngot:\n%v\nwant:\n%v", mat.Formatted(&got), mat.Formatted(&want)) } } } func TestMatCloneFrom(t *testing.T) { + const tol = 1e-16 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") + want := randomMat(rnd) + got := NewMat(nil) + got.CloneFrom(want) + if !mat.EqualApprox(got, want, tol) { + 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 got != expect { - t.Error("r3.Cross(v1,v2) not match with r3.Skew(v1)*v2") + if d := want.Sub(got); d.Dot(d) > tol { + t.Errorf("r3.Cross(v1,v2) does not agree with r3.Skew(v1)*v2: got:%v want:%v", got, want) } } } 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 { - t.Error("VecDense.MulVec(dense.T()) not equal to r3.Mat.MulVec(r3.Vec)") + want := Vec{X: vd.AtVec(0), Y: vd.AtVec(1), Z: vd.AtVec(2)} + 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 diff --git a/spatial/r3/mat_unsafe.go b/spatial/r3/mat_unsafe.go new file mode 100644 index 00000000..2281649d --- /dev/null +++ b/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))[:] +}