--- /dev/null
+// 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"
+)
+
+var (
+ dense *Dense
+
+ _ Matrix = dense
+ _ Mutable = dense
+
+ _ Cloner = dense
+ _ RowViewer = dense
+ _ ColViewer = dense
+ _ RawRowViewer = dense
+ _ Grower = dense
+
+ _ RawMatrixSetter = dense
+ _ RawMatrixer = dense
+
+ _ Reseter = dense
+)
+
+// Dense is a dense matrix representation.
+type Dense struct {
+ mat blas64.General
+
+ capRows, capCols int
+}
+
+// NewDense creates a new Dense matrix with r rows and c columns. If data == nil,
+// a new slice is allocated for the backing slice. If len(data) == r*c, data is
+// used as the backing slice, and changes to the elements of the returned Dense
+// will be reflected in data. If neither of these is true, NewDense will panic.
+//
+// The data must be arranged in row-major order, i.e. the (i*c + j)-th
+// element in the data slice is the {i, j}-th element in the matrix.
+func NewDense(r, c int, data []float64) *Dense {
+ if data != nil && r*c != len(data) {
+ panic(ErrShape)
+ }
+ if data == nil {
+ data = make([]float64, r*c)
+ }
+ return &Dense{
+ mat: blas64.General{
+ Rows: r,
+ Cols: c,
+ Stride: c,
+ Data: data,
+ },
+ capRows: r,
+ capCols: c,
+ }
+}
+
+// reuseAs resizes an empty matrix to a r×c matrix,
+// or checks that a non-empty matrix is r×c.
+//
+// reuseAs must be kept in sync with reuseAsZeroed.
+func (m *Dense) reuseAs(r, c int) {
+ if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
+ // Panic as a string, not a mat.Error.
+ panic("mat: caps not correctly set")
+ }
+ if m.IsZero() {
+ m.mat = blas64.General{
+ Rows: r,
+ Cols: c,
+ Stride: c,
+ Data: use(m.mat.Data, r*c),
+ }
+ m.capRows = r
+ m.capCols = c
+ return
+ }
+ if r != m.mat.Rows || c != m.mat.Cols {
+ panic(ErrShape)
+ }
+}
+
+// reuseAsZeroed resizes an empty matrix to a r×c matrix,
+// or checks that a non-empty matrix is r×c. It zeroes
+// all the elements of the matrix.
+//
+// reuseAsZeroed must be kept in sync with reuseAs.
+func (m *Dense) reuseAsZeroed(r, c int) {
+ if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
+ // Panic as a string, not a mat.Error.
+ panic("mat: caps not correctly set")
+ }
+ if m.IsZero() {
+ m.mat = blas64.General{
+ Rows: r,
+ Cols: c,
+ Stride: c,
+ Data: useZeroed(m.mat.Data, r*c),
+ }
+ m.capRows = r
+ m.capCols = c
+ return
+ }
+ if r != m.mat.Rows || c != m.mat.Cols {
+ panic(ErrShape)
+ }
+ for i := 0; i < r; i++ {
+ zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c])
+ }
+}
+
+// untranspose untransposes a matrix if applicable. If a is an Untransposer, then
+// untranspose returns the underlying matrix and true. If it is not, then it returns
+// the input matrix and false.
+func untranspose(a Matrix) (Matrix, bool) {
+ if ut, ok := a.(Untransposer); ok {
+ return ut.Untranspose(), true
+ }
+ return a, false
+}
+
+// isolatedWorkspace returns a new dense matrix w with the size of a and
+// returns a callback to defer which performs cleanup at the return of the call.
+// This should be used when a method receiver is the same pointer as an input argument.
+func (m *Dense) isolatedWorkspace(a Matrix) (w *Dense, restore func()) {
+ r, c := a.Dims()
+ w = getWorkspace(r, c, false)
+ return w, func() {
+ m.Copy(w)
+ putWorkspace(w)
+ }
+}
+
+// Reset zeros the dimensions of the matrix so that it can be reused as the
+// receiver of a dimensionally restricted operation.
+//
+// See the Reseter interface for more information.
+func (m *Dense) Reset() {
+ // Row, Cols and Stride must be zeroed in unison.
+ m.mat.Rows, m.mat.Cols, m.mat.Stride = 0, 0, 0
+ m.capRows, m.capCols = 0, 0
+ m.mat.Data = m.mat.Data[:0]
+}
+
+// IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
+// receiver for size-restricted operations. Dense matrices can be zeroed using Reset.
+func (m *Dense) IsZero() bool {
+ // It must be the case that m.Dims() returns
+ // zeros in this case. See comment in Reset().
+ return m.mat.Stride == 0
+}
+
+// asTriDense returns a TriDense with the given size and side. The backing data
+// of the TriDense is the same as the receiver.
+func (m *Dense) asTriDense(n int, diag blas.Diag, uplo blas.Uplo) *TriDense {
+ return &TriDense{
+ mat: blas64.Triangular{
+ N: n,
+ Stride: m.mat.Stride,
+ Data: m.mat.Data,
+ Uplo: uplo,
+ Diag: diag,
+ },
+ cap: n,
+ }
+}
+
+// DenseCopyOf returns a newly allocated copy of the elements of a.
+func DenseCopyOf(a Matrix) *Dense {
+ d := &Dense{}
+ d.Clone(a)
+ return d
+}
+
+// SetRawMatrix sets the underlying blas64.General used by the receiver.
+// Changes to elements in the receiver following the call will be reflected
+// in b.
+func (m *Dense) SetRawMatrix(b blas64.General) {
+ m.capRows, m.capCols = b.Rows, b.Cols
+ m.mat = b
+}
+
+// RawMatrix returns the underlying blas64.General used by the receiver.
+// Changes to elements in the receiver following the call will be reflected
+// in returned blas64.General.
+func (m *Dense) RawMatrix() blas64.General { return m.mat }
+
+// Dims returns the number of rows and columns in the matrix.
+func (m *Dense) Dims() (r, c int) { return m.mat.Rows, m.mat.Cols }
+
+// Caps returns the number of rows and columns in the backing matrix.
+func (m *Dense) Caps() (r, c int) { return m.capRows, m.capCols }
+
+// T performs an implicit transpose by returning the receiver inside a Transpose.
+func (m *Dense) T() Matrix {
+ return Transpose{m}
+}
+
+// ColView returns a Vector reflecting the column j, backed by the matrix data.
+//
+// See ColViewer for more information.
+func (m *Dense) ColView(j int) Vector {
+ var v VecDense
+ v.ColViewOf(m, j)
+ return &v
+}
+
+// SetCol sets the values in the specified column of the matrix to the values
+// in src. len(src) must equal the number of rows in the receiver.
+func (m *Dense) SetCol(j int, src []float64) {
+ if j >= m.mat.Cols || j < 0 {
+ panic(ErrColAccess)
+ }
+ if len(src) != m.mat.Rows {
+ panic(ErrColLength)
+ }
+
+ blas64.Copy(m.mat.Rows,
+ blas64.Vector{Inc: 1, Data: src},
+ blas64.Vector{Inc: m.mat.Stride, Data: m.mat.Data[j:]},
+ )
+}
+
+// SetRow sets the values in the specified rows of the matrix to the values
+// in src. len(src) must equal the number of columns in the receiver.
+func (m *Dense) SetRow(i int, src []float64) {
+ if i >= m.mat.Rows || i < 0 {
+ panic(ErrRowAccess)
+ }
+ if len(src) != m.mat.Cols {
+ panic(ErrRowLength)
+ }
+
+ copy(m.rawRowView(i), src)
+}
+
+// RowView returns row i of the matrix data represented as a column vector,
+// backed by the matrix data.
+//
+// See RowViewer for more information.
+func (m *Dense) RowView(i int) Vector {
+ var v VecDense
+ v.RowViewOf(m, i)
+ return &v
+}
+
+// RawRowView returns a slice backed by the same array as backing the
+// receiver.
+func (m *Dense) RawRowView(i int) []float64 {
+ if i >= m.mat.Rows || i < 0 {
+ panic(ErrRowAccess)
+ }
+ return m.rawRowView(i)
+}
+
+func (m *Dense) rawRowView(i int) []float64 {
+ return m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+m.mat.Cols]
+}
+
+// Slice returns a new Matrix that shares backing data with the receiver.
+// The returned matrix starts at {i,j} of the receiver and extends k-i rows
+// and l-j columns. The final row in the resulting matrix is k-1 and the
+// final column is l-1.
+// Slice panics with ErrIndexOutOfRange if the slice is outside the capacity
+// of the receiver.
+func (m *Dense) Slice(i, k, j, l int) Matrix {
+ mr, mc := m.Caps()
+ if i < 0 || mr <= i || j < 0 || mc <= j || k <= i || mr < k || l <= j || mc < l {
+ panic(ErrIndexOutOfRange)
+ }
+ t := *m
+ t.mat.Data = t.mat.Data[i*t.mat.Stride+j : (k-1)*t.mat.Stride+l]
+ t.mat.Rows = k - i
+ t.mat.Cols = l - j
+ t.capRows -= i
+ t.capCols -= j
+ return &t
+}
+
+// Grow returns the receiver expanded by r rows and c columns. If the dimensions
+// of the expanded matrix are outside the capacities of the receiver a new
+// allocation is made, otherwise not. Note the receiver itself is not modified
+// during the call to Grow.
+func (m *Dense) Grow(r, c int) Matrix {
+ if r < 0 || c < 0 {
+ panic(ErrIndexOutOfRange)
+ }
+ if r == 0 && c == 0 {
+ return m
+ }
+
+ r += m.mat.Rows
+ c += m.mat.Cols
+
+ var t Dense
+ switch {
+ case m.mat.Rows == 0 || m.mat.Cols == 0:
+ t.mat = blas64.General{
+ Rows: r,
+ Cols: c,
+ Stride: c,
+ // We zero because we don't know how the matrix will be used.
+ // In other places, the mat is immediately filled with a result;
+ // this is not the case here.
+ Data: useZeroed(m.mat.Data, r*c),
+ }
+ case r > m.capRows || c > m.capCols:
+ cr := max(r, m.capRows)
+ cc := max(c, m.capCols)
+ t.mat = blas64.General{
+ Rows: r,
+ Cols: c,
+ Stride: cc,
+ Data: make([]float64, cr*cc),
+ }
+ t.capRows = cr
+ t.capCols = cc
+ // Copy the complete matrix over to the new matrix.
+ // Including elements not currently visible. Use a temporary structure
+ // to avoid modifying the receiver.
+ var tmp Dense
+ tmp.mat = blas64.General{
+ Rows: m.mat.Rows,
+ Cols: m.mat.Cols,
+ Stride: m.mat.Stride,
+ Data: m.mat.Data,
+ }
+ tmp.capRows = m.capRows
+ tmp.capCols = m.capCols
+ t.Copy(&tmp)
+ return &t
+ default:
+ t.mat = blas64.General{
+ Data: m.mat.Data[:(r-1)*m.mat.Stride+c],
+ Rows: r,
+ Cols: c,
+ Stride: m.mat.Stride,
+ }
+ }
+ t.capRows = r
+ t.capCols = c
+ return &t
+}
+
+// Clone makes a copy of a into the receiver, overwriting the previous value of
+// the receiver. The clone operation does not make any restriction on shape and
+// will not cause shadowing.
+//
+// See the Cloner interface for more information.
+func (m *Dense) Clone(a Matrix) {
+ r, c := a.Dims()
+ mat := blas64.General{
+ Rows: r,
+ Cols: c,
+ Stride: c,
+ }
+ m.capRows, m.capCols = r, c
+
+ aU, trans := untranspose(a)
+ switch aU := aU.(type) {
+ case RawMatrixer:
+ amat := aU.RawMatrix()
+ mat.Data = make([]float64, r*c)
+ if trans {
+ for i := 0; i < r; i++ {
+ blas64.Copy(c,
+ blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
+ blas64.Vector{Inc: 1, Data: mat.Data[i*c : (i+1)*c]})
+ }
+ } else {
+ for i := 0; i < r; i++ {
+ copy(mat.Data[i*c:(i+1)*c], amat.Data[i*amat.Stride:i*amat.Stride+c])
+ }
+ }
+ case *VecDense:
+ amat := aU.mat
+ mat.Data = make([]float64, aU.n)
+ blas64.Copy(aU.n,
+ blas64.Vector{Inc: amat.Inc, Data: amat.Data},
+ blas64.Vector{Inc: 1, Data: mat.Data})
+ default:
+ mat.Data = make([]float64, r*c)
+ w := *m
+ w.mat = mat
+ for i := 0; i < r; i++ {
+ for j := 0; j < c; j++ {
+ w.set(i, j, a.At(i, j))
+ }
+ }
+ *m = w
+ return
+ }
+ m.mat = mat
+}
+
+// Copy 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 matrices and
+// returns the number of rows and columns it copied. If a aliases the receiver
+// and is a transposed Dense or VecDense, with a non-unitary increment, Copy will
+// panic.
+//
+// See the Copier interface for more information.
+func (m *Dense) Copy(a Matrix) (r, c int) {
+ r, c = a.Dims()
+ if a == m {
+ return r, c
+ }
+ r = min(r, m.mat.Rows)
+ c = min(c, m.mat.Cols)
+ if r == 0 || c == 0 {
+ return 0, 0
+ }
+
+ aU, trans := untranspose(a)
+ switch aU := aU.(type) {
+ case RawMatrixer:
+ amat := aU.RawMatrix()
+ if trans {
+ if amat.Stride != 1 {
+ m.checkOverlap(amat)
+ }
+ for i := 0; i < r; i++ {
+ blas64.Copy(c,
+ blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
+ blas64.Vector{Inc: 1, Data: m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c]})
+ }
+ } else {
+ switch o := offset(m.mat.Data, amat.Data); {
+ case o < 0:
+ for i := r - 1; i >= 0; i-- {
+ copy(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], amat.Data[i*amat.Stride:i*amat.Stride+c])
+ }
+ case o > 0:
+ for i := 0; i < r; i++ {
+ copy(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], amat.Data[i*amat.Stride:i*amat.Stride+c])
+ }
+ default:
+ // Nothing to do.
+ }
+ }
+ case *VecDense:
+ var n, stride int
+ amat := aU.mat
+ if trans {
+ if amat.Inc != 1 {
+ m.checkOverlap(aU.asGeneral())
+ }
+ n = c
+ stride = 1
+ } else {
+ n = r
+ stride = m.mat.Stride
+ }
+ if amat.Inc == 1 && stride == 1 {
+ copy(m.mat.Data, amat.Data[:n])
+ break
+ }
+ switch o := offset(m.mat.Data, amat.Data); {
+ case o < 0:
+ blas64.Copy(n,
+ blas64.Vector{Inc: -amat.Inc, Data: amat.Data},
+ blas64.Vector{Inc: -stride, Data: m.mat.Data})
+ case o > 0:
+ blas64.Copy(n,
+ blas64.Vector{Inc: amat.Inc, Data: amat.Data},
+ blas64.Vector{Inc: stride, Data: m.mat.Data})
+ default:
+ // Nothing to do.
+ }
+ default:
+ for i := 0; i < r; i++ {
+ for j := 0; j < c; j++ {
+ m.set(i, j, a.At(i, j))
+ }
+ }
+ }
+
+ return r, c
+}
+
+// Stack appends the rows of b onto the rows of a, placing the result into the
+// receiver with b placed in the greater indexed rows. Stack will panic if the
+// two input matrices do not have the same number of columns or the constructed
+// stacked matrix is not the same shape as the receiver.
+func (m *Dense) Stack(a, b Matrix) {
+ ar, ac := a.Dims()
+ br, bc := b.Dims()
+ if ac != bc || m == a || m == b {
+ panic(ErrShape)
+ }
+
+ m.reuseAs(ar+br, ac)
+
+ m.Copy(a)
+ w := m.Slice(ar, ar+br, 0, bc).(*Dense)
+ w.Copy(b)
+}
+
+// Augment creates the augmented matrix of a and b, where b is placed in the
+// greater indexed columns. Augment will panic if the two input matrices do
+// not have the same number of rows or the constructed augmented matrix is
+// not the same shape as the receiver.
+func (m *Dense) Augment(a, b Matrix) {
+ ar, ac := a.Dims()
+ br, bc := b.Dims()
+ if ar != br || m == a || m == b {
+ panic(ErrShape)
+ }
+
+ m.reuseAs(ar, ac+bc)
+
+ m.Copy(a)
+ w := m.Slice(0, br, ac, ac+bc).(*Dense)
+ w.Copy(b)
+}
OSDN Git Service
