k3s/vendor/gonum.org/v1/gonum/mat/vector.go

742 lines
16 KiB
Go
Raw Blame History

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

// Copyright ©2013 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 mat
import (
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
"gonum.org/v1/gonum/internal/asm/f64"
)
var (
vector *VecDense
_ Matrix = vector
_ Vector = vector
_ Reseter = vector
)
// Vector is a vector.
type Vector interface {
Matrix
AtVec(int) float64
Len() int
}
// TransposeVec is a type for performing an implicit transpose of a Vector.
// It implements the Vector interface, returning values from the transpose
// of the vector within.
type TransposeVec struct {
Vector Vector
}
// At returns the value of the element at row i and column j of the transposed
// matrix, that is, row j and column i of the Vector field.
func (t TransposeVec) At(i, j int) float64 {
return t.Vector.At(j, i)
}
// AtVec returns the element at position i. It panics if i is out of bounds.
func (t TransposeVec) AtVec(i int) float64 {
return t.Vector.AtVec(i)
}
// Dims returns the dimensions of the transposed vector.
func (t TransposeVec) Dims() (r, c int) {
c, r = t.Vector.Dims()
return r, c
}
// T performs an implicit transpose by returning the Vector field.
func (t TransposeVec) T() Matrix {
return t.Vector
}
// Len returns the number of columns in the vector.
func (t TransposeVec) Len() int {
return t.Vector.Len()
}
// TVec performs an implicit transpose by returning the Vector field.
func (t TransposeVec) TVec() Vector {
return t.Vector
}
// Untranspose returns the Vector field.
func (t TransposeVec) Untranspose() Matrix {
return t.Vector
}
func (t TransposeVec) UntransposeVec() Vector {
return t.Vector
}
// VecDense represents a column vector.
type VecDense struct {
mat blas64.Vector
// A BLAS vector can have a negative increment, but allowing this
// in the mat type complicates a lot of code, and doesn't gain anything.
// VecDense must have positive increment in this package.
}
// NewVecDense creates a new VecDense of length n. If data == nil,
// a new slice is allocated for the backing slice. If len(data) == n, data is
// used as the backing slice, and changes to the elements of the returned VecDense
// will be reflected in data. If neither of these is true, NewVecDense will panic.
// NewVecDense will panic if n is zero.
func NewVecDense(n int, data []float64) *VecDense {
if n <= 0 {
if n == 0 {
panic(ErrZeroLength)
}
panic("mat: negative dimension")
}
if len(data) != n && data != nil {
panic(ErrShape)
}
if data == nil {
data = make([]float64, n)
}
return &VecDense{
mat: blas64.Vector{
N: n,
Inc: 1,
Data: data,
},
}
}
// SliceVec returns a new Vector that shares backing data with the receiver.
// The returned matrix starts at i of the receiver and extends k-i elements.
// SliceVec panics with ErrIndexOutOfRange if the slice is outside the capacity
// of the receiver.
func (v *VecDense) SliceVec(i, k int) Vector {
if i < 0 || k <= i || v.Cap() < k {
panic(ErrIndexOutOfRange)
}
return &VecDense{
mat: blas64.Vector{
N: k - i,
Inc: v.mat.Inc,
Data: v.mat.Data[i*v.mat.Inc : (k-1)*v.mat.Inc+1],
},
}
}
// Dims returns the number of rows and columns in the matrix. Columns is always 1
// for a non-Reset vector.
func (v *VecDense) Dims() (r, c int) {
if v.IsZero() {
return 0, 0
}
return v.mat.N, 1
}
// Caps returns the number of rows and columns in the backing matrix. Columns is always 1
// for a non-Reset vector.
func (v *VecDense) Caps() (r, c int) {
if v.IsZero() {
return 0, 0
}
return v.Cap(), 1
}
// Len returns the length of the vector.
func (v *VecDense) Len() int {
return v.mat.N
}
// Cap returns the capacity of the vector.
func (v *VecDense) Cap() int {
if v.IsZero() {
return 0
}
return (cap(v.mat.Data)-1)/v.mat.Inc + 1
}
// T performs an implicit transpose by returning the receiver inside a Transpose.
func (v *VecDense) T() Matrix {
return Transpose{v}
}
// TVec performs an implicit transpose by returning the receiver inside a TransposeVec.
func (v *VecDense) TVec() Vector {
return TransposeVec{v}
}
// Reset zeros the length of the vector so that it can be reused as the
// receiver of a dimensionally restricted operation.
//
// See the Reseter interface for more information.
func (v *VecDense) Reset() {
// No change of Inc or N to 0 may be
// made unless both are set to 0.
v.mat.Inc = 0
v.mat.N = 0
v.mat.Data = v.mat.Data[:0]
}
// Zero sets all of the matrix elements to zero.
func (v *VecDense) Zero() {
for i := 0; i < v.mat.N; i++ {
v.mat.Data[v.mat.Inc*i] = 0
}
}
// CloneVec makes a copy of a into the receiver, overwriting the previous value
// of the receiver.
func (v *VecDense) CloneVec(a Vector) {
if v == a {
return
}
n := a.Len()
v.mat = blas64.Vector{
N: n,
Inc: 1,
Data: use(v.mat.Data, n),
}
if r, ok := a.(RawVectorer); ok {
blas64.Copy(r.RawVector(), v.mat)
return
}
for i := 0; i < a.Len(); i++ {
v.SetVec(i, a.AtVec(i))
}
}
// VecDenseCopyOf returns a newly allocated copy of the elements of a.
func VecDenseCopyOf(a Vector) *VecDense {
v := &VecDense{}
v.CloneVec(a)
return v
}
func (v *VecDense) RawVector() blas64.Vector {
return v.mat
}
// CopyVec makes a copy of elements of a into the receiver. It is similar to the
// built-in copy; it copies as much as the overlap between the two vectors and
// returns the number of elements it copied.
func (v *VecDense) CopyVec(a Vector) int {
n := min(v.Len(), a.Len())
if v == a {
return n
}
if r, ok := a.(RawVectorer); ok {
blas64.Copy(r.RawVector(), v.mat)
return n
}
for i := 0; i < n; i++ {
v.setVec(i, a.AtVec(i))
}
return n
}
// ScaleVec scales the vector a by alpha, placing the result in the receiver.
func (v *VecDense) ScaleVec(alpha float64, a Vector) {
n := a.Len()
if v == a {
if v.mat.Inc == 1 {
f64.ScalUnitary(alpha, v.mat.Data)
return
}
f64.ScalInc(alpha, v.mat.Data, uintptr(n), uintptr(v.mat.Inc))
return
}
v.reuseAs(n)
if rv, ok := a.(RawVectorer); ok {
mat := rv.RawVector()
v.checkOverlap(mat)
if v.mat.Inc == 1 && mat.Inc == 1 {
f64.ScalUnitaryTo(v.mat.Data, alpha, mat.Data)
return
}
f64.ScalIncTo(v.mat.Data, uintptr(v.mat.Inc),
alpha, mat.Data, uintptr(n), uintptr(mat.Inc))
return
}
for i := 0; i < n; i++ {
v.setVec(i, alpha*a.AtVec(i))
}
}
// AddScaledVec adds the vectors a and alpha*b, placing the result in the receiver.
func (v *VecDense) AddScaledVec(a Vector, alpha float64, b Vector) {
if alpha == 1 {
v.AddVec(a, b)
return
}
if alpha == -1 {
v.SubVec(a, b)
return
}
ar := a.Len()
br := b.Len()
if ar != br {
panic(ErrShape)
}
var amat, bmat blas64.Vector
fast := true
aU, _ := untranspose(a)
if rv, ok := aU.(RawVectorer); ok {
amat = rv.RawVector()
if v != a {
v.checkOverlap(amat)
}
} else {
fast = false
}
bU, _ := untranspose(b)
if rv, ok := bU.(RawVectorer); ok {
bmat = rv.RawVector()
if v != b {
v.checkOverlap(bmat)
}
} else {
fast = false
}
v.reuseAs(ar)
switch {
case alpha == 0: // v <- a
if v == a {
return
}
v.CopyVec(a)
case v == a && v == b: // v <- v + alpha * v = (alpha + 1) * v
blas64.Scal(alpha+1, v.mat)
case !fast: // v <- a + alpha * b without blas64 support.
for i := 0; i < ar; i++ {
v.setVec(i, a.AtVec(i)+alpha*b.AtVec(i))
}
case v == a && v != b: // v <- v + alpha * b
if v.mat.Inc == 1 && bmat.Inc == 1 {
// Fast path for a common case.
f64.AxpyUnitaryTo(v.mat.Data, alpha, bmat.Data, amat.Data)
} else {
f64.AxpyInc(alpha, bmat.Data, v.mat.Data,
uintptr(ar), uintptr(bmat.Inc), uintptr(v.mat.Inc), 0, 0)
}
default: // v <- a + alpha * b or v <- a + alpha * v
if v.mat.Inc == 1 && amat.Inc == 1 && bmat.Inc == 1 {
// Fast path for a common case.
f64.AxpyUnitaryTo(v.mat.Data, alpha, bmat.Data, amat.Data)
} else {
f64.AxpyIncTo(v.mat.Data, uintptr(v.mat.Inc), 0,
alpha, bmat.Data, amat.Data,
uintptr(ar), uintptr(bmat.Inc), uintptr(amat.Inc), 0, 0)
}
}
}
// AddVec adds the vectors a and b, placing the result in the receiver.
func (v *VecDense) AddVec(a, b Vector) {
ar := a.Len()
br := b.Len()
if ar != br {
panic(ErrShape)
}
v.reuseAs(ar)
aU, _ := untranspose(a)
bU, _ := untranspose(b)
if arv, ok := aU.(RawVectorer); ok {
if brv, ok := bU.(RawVectorer); ok {
amat := arv.RawVector()
bmat := brv.RawVector()
if v != a {
v.checkOverlap(amat)
}
if v != b {
v.checkOverlap(bmat)
}
if v.mat.Inc == 1 && amat.Inc == 1 && bmat.Inc == 1 {
// Fast path for a common case.
f64.AxpyUnitaryTo(v.mat.Data, 1, bmat.Data, amat.Data)
return
}
f64.AxpyIncTo(v.mat.Data, uintptr(v.mat.Inc), 0,
1, bmat.Data, amat.Data,
uintptr(ar), uintptr(bmat.Inc), uintptr(amat.Inc), 0, 0)
return
}
}
for i := 0; i < ar; i++ {
v.setVec(i, a.AtVec(i)+b.AtVec(i))
}
}
// SubVec subtracts the vector b from a, placing the result in the receiver.
func (v *VecDense) SubVec(a, b Vector) {
ar := a.Len()
br := b.Len()
if ar != br {
panic(ErrShape)
}
v.reuseAs(ar)
aU, _ := untranspose(a)
bU, _ := untranspose(b)
if arv, ok := aU.(RawVectorer); ok {
if brv, ok := bU.(RawVectorer); ok {
amat := arv.RawVector()
bmat := brv.RawVector()
if v != a {
v.checkOverlap(amat)
}
if v != b {
v.checkOverlap(bmat)
}
if v.mat.Inc == 1 && amat.Inc == 1 && bmat.Inc == 1 {
// Fast path for a common case.
f64.AxpyUnitaryTo(v.mat.Data, -1, bmat.Data, amat.Data)
return
}
f64.AxpyIncTo(v.mat.Data, uintptr(v.mat.Inc), 0,
-1, bmat.Data, amat.Data,
uintptr(ar), uintptr(bmat.Inc), uintptr(amat.Inc), 0, 0)
return
}
}
for i := 0; i < ar; i++ {
v.setVec(i, a.AtVec(i)-b.AtVec(i))
}
}
// MulElemVec performs element-wise multiplication of a and b, placing the result
// in the receiver.
func (v *VecDense) MulElemVec(a, b Vector) {
ar := a.Len()
br := b.Len()
if ar != br {
panic(ErrShape)
}
v.reuseAs(ar)
aU, _ := untranspose(a)
bU, _ := untranspose(b)
if arv, ok := aU.(RawVectorer); ok {
if brv, ok := bU.(RawVectorer); ok {
amat := arv.RawVector()
bmat := brv.RawVector()
if v != a {
v.checkOverlap(amat)
}
if v != b {
v.checkOverlap(bmat)
}
if v.mat.Inc == 1 && amat.Inc == 1 && bmat.Inc == 1 {
// Fast path for a common case.
for i, a := range amat.Data {
v.mat.Data[i] = a * bmat.Data[i]
}
return
}
var ia, ib int
for i := 0; i < ar; i++ {
v.setVec(i, amat.Data[ia]*bmat.Data[ib])
ia += amat.Inc
ib += bmat.Inc
}
return
}
}
for i := 0; i < ar; i++ {
v.setVec(i, a.AtVec(i)*b.AtVec(i))
}
}
// DivElemVec performs element-wise division of a by b, placing the result
// in the receiver.
func (v *VecDense) DivElemVec(a, b Vector) {
ar := a.Len()
br := b.Len()
if ar != br {
panic(ErrShape)
}
v.reuseAs(ar)
aU, _ := untranspose(a)
bU, _ := untranspose(b)
if arv, ok := aU.(RawVectorer); ok {
if brv, ok := bU.(RawVectorer); ok {
amat := arv.RawVector()
bmat := brv.RawVector()
if v != a {
v.checkOverlap(amat)
}
if v != b {
v.checkOverlap(bmat)
}
if v.mat.Inc == 1 && amat.Inc == 1 && bmat.Inc == 1 {
// Fast path for a common case.
for i, a := range amat.Data {
v.setVec(i, a/bmat.Data[i])
}
return
}
var ia, ib int
for i := 0; i < ar; i++ {
v.setVec(i, amat.Data[ia]/bmat.Data[ib])
ia += amat.Inc
ib += bmat.Inc
}
}
}
for i := 0; i < ar; i++ {
v.setVec(i, a.AtVec(i)/b.AtVec(i))
}
}
// MulVec computes a * b. The result is stored into the receiver.
// MulVec panics if the number of columns in a does not equal the number of rows in b
// or if the number of columns in b does not equal 1.
func (v *VecDense) MulVec(a Matrix, b Vector) {
r, c := a.Dims()
br, bc := b.Dims()
if c != br || bc != 1 {
panic(ErrShape)
}
aU, trans := untranspose(a)
var bmat blas64.Vector
fast := true
bU, _ := untranspose(b)
if rv, ok := bU.(RawVectorer); ok {
bmat = rv.RawVector()
if v != b {
v.checkOverlap(bmat)
}
} else {
fast = false
}
v.reuseAs(r)
var restore func()
if v == aU {
v, restore = v.isolatedWorkspace(aU.(*VecDense))
defer restore()
} else if v == b {
v, restore = v.isolatedWorkspace(b)
defer restore()
}
// TODO(kortschak): Improve the non-fast paths.
switch aU := aU.(type) {
case Vector:
if b.Len() == 1 {
// {n,1} x {1,1}
v.ScaleVec(b.AtVec(0), aU)
return
}
// {1,n} x {n,1}
if fast {
if rv, ok := aU.(RawVectorer); ok {
amat := rv.RawVector()
if v != aU {
v.checkOverlap(amat)
}
if amat.Inc == 1 && bmat.Inc == 1 {
// Fast path for a common case.
v.setVec(0, f64.DotUnitary(amat.Data, bmat.Data))
return
}
v.setVec(0, f64.DotInc(amat.Data, bmat.Data,
uintptr(c), uintptr(amat.Inc), uintptr(bmat.Inc), 0, 0))
return
}
}
var sum float64
for i := 0; i < c; i++ {
sum += aU.AtVec(i) * b.AtVec(i)
}
v.setVec(0, sum)
return
case RawSymmetricer:
if fast {
amat := aU.RawSymmetric()
// We don't know that a is a *SymDense, so make
// a temporary SymDense to check overlap.
(&SymDense{mat: amat}).checkOverlap(v.asGeneral())
blas64.Symv(1, amat, bmat, 0, v.mat)
return
}
case RawTriangular:
v.CopyVec(b)
amat := aU.RawTriangular()
// We don't know that a is a *TriDense, so make
// a temporary TriDense to check overlap.
(&TriDense{mat: amat}).checkOverlap(v.asGeneral())
ta := blas.NoTrans
if trans {
ta = blas.Trans
}
blas64.Trmv(ta, amat, v.mat)
case RawMatrixer:
if fast {
amat := aU.RawMatrix()
// We don't know that a is a *Dense, so make
// a temporary Dense to check overlap.
(&Dense{mat: amat}).checkOverlap(v.asGeneral())
t := blas.NoTrans
if trans {
t = blas.Trans
}
blas64.Gemv(t, 1, amat, bmat, 0, v.mat)
return
}
default:
if fast {
for i := 0; i < r; i++ {
var f float64
for j := 0; j < c; j++ {
f += a.At(i, j) * bmat.Data[j*bmat.Inc]
}
v.setVec(i, f)
}
return
}
}
for i := 0; i < r; i++ {
var f float64
for j := 0; j < c; j++ {
f += a.At(i, j) * b.AtVec(j)
}
v.setVec(i, f)
}
}
// reuseAs resizes an empty vector to a r×1 vector,
// or checks that a non-empty matrix is r×1.
func (v *VecDense) reuseAs(r int) {
if r == 0 {
panic(ErrZeroLength)
}
if v.IsZero() {
v.mat = blas64.Vector{
N: r,
Inc: 1,
Data: use(v.mat.Data, r),
}
return
}
if r != v.mat.N {
panic(ErrShape)
}
}
// IsZero returns whether the receiver is zero-sized. Zero-sized vectors can be the
// receiver for size-restricted operations. VecDenses can be zeroed using Reset.
func (v *VecDense) IsZero() bool {
// It must be the case that v.Dims() returns
// zeros in this case. See comment in Reset().
return v.mat.Inc == 0
}
func (v *VecDense) isolatedWorkspace(a Vector) (n *VecDense, restore func()) {
l := a.Len()
if l == 0 {
panic(ErrZeroLength)
}
n = getWorkspaceVec(l, false)
return n, func() {
v.CopyVec(n)
putWorkspaceVec(n)
}
}
// asDense returns a Dense representation of the receiver with the same
// underlying data.
func (v *VecDense) asDense() *Dense {
return &Dense{
mat: v.asGeneral(),
capRows: v.mat.N,
capCols: 1,
}
}
// asGeneral returns a blas64.General representation of the receiver with the
// same underlying data.
func (v *VecDense) asGeneral() blas64.General {
return blas64.General{
Rows: v.mat.N,
Cols: 1,
Stride: v.mat.Inc,
Data: v.mat.Data,
}
}
// ColViewOf reflects the column j of the RawMatrixer m, into the receiver
// backed by the same underlying data. The length of the receiver must either be
// zero or match the number of rows in m.
func (v *VecDense) ColViewOf(m RawMatrixer, j int) {
rm := m.RawMatrix()
if j >= rm.Cols || j < 0 {
panic(ErrColAccess)
}
if !v.IsZero() && v.mat.N != rm.Rows {
panic(ErrShape)
}
v.mat.Inc = rm.Stride
v.mat.Data = rm.Data[j : (rm.Rows-1)*rm.Stride+j+1]
v.mat.N = rm.Rows
}
// RowViewOf reflects the row i of the RawMatrixer m, into the receiver
// backed by the same underlying data. The length of the receiver must either be
// zero or match the number of columns in m.
func (v *VecDense) RowViewOf(m RawMatrixer, i int) {
rm := m.RawMatrix()
if i >= rm.Rows || i < 0 {
panic(ErrRowAccess)
}
if !v.IsZero() && v.mat.N != rm.Cols {
panic(ErrShape)
}
v.mat.Inc = 1
v.mat.Data = rm.Data[i*rm.Stride : i*rm.Stride+rm.Cols]
v.mat.N = rm.Cols
}