vendor/gonum.org/v1/gonum/mat/dense.go

   1 // Copyright ©2013 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package mat
   6
   7 import (
   8         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/blas/blas64"
  10 )
  11
  12 var (
  13         dense *Dense
  14
  15         _ Matrix  = dense
  16         _ Mutable = dense
  17
  18         _ Cloner       = dense
  19         _ RowViewer    = dense
  20         _ ColViewer    = dense
  21         _ RawRowViewer = dense
  22         _ Grower       = dense
  23
  24         _ RawMatrixSetter = dense
  25         _ RawMatrixer     = dense
  26
  27         _ Reseter = dense
  28 )
  29
  30 // Dense is a dense matrix representation.
  31 type Dense struct {
  32         mat blas64.General
  33
  34         capRows, capCols int
  35 }
  36
  37 // NewDense creates a new Dense matrix with r rows and c columns. If data == nil,
  38 // a new slice is allocated for the backing slice. If len(data) == r*c, data is
  39 // used as the backing slice, and changes to the elements of the returned Dense
  40 // will be reflected in data. If neither of these is true, NewDense will panic.
  41 //
  42 // The data must be arranged in row-major order, i.e. the (i*c + j)-th
  43 // element in the data slice is the {i, j}-th element in the matrix.
  44 func NewDense(r, c int, data []float64) *Dense {
  45         if data != nil && r*c != len(data) {
  46                 panic(ErrShape)
  47         }
  48         if data == nil {
  49                 data = make([]float64, r*c)
  50         }
  51         return &Dense{
  52                 mat: blas64.General{
  53                         Rows:   r,
  54                         Cols:   c,
  55                         Stride: c,
  56                         Data:   data,
  57                 },
  58                 capRows: r,
  59                 capCols: c,
  60         }
  61 }
  62
  63 // reuseAs resizes an empty matrix to a r×c matrix,
  64 // or checks that a non-empty matrix is r×c.
  65 //
  66 // reuseAs must be kept in sync with reuseAsZeroed.
  67 func (m *Dense) reuseAs(r, c int) {
  68         if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
  69                 // Panic as a string, not a mat.Error.
  70                 panic("mat: caps not correctly set")
  71         }
  72         if m.IsZero() {
  73                 m.mat = blas64.General{
  74                         Rows:   r,
  75                         Cols:   c,
  76                         Stride: c,
  77                         Data:   use(m.mat.Data, r*c),
  78                 }
  79                 m.capRows = r
  80                 m.capCols = c
  81                 return
  82         }
  83         if r != m.mat.Rows || c != m.mat.Cols {
  84                 panic(ErrShape)
  85         }
  86 }
  87
  88 // reuseAsZeroed resizes an empty matrix to a r×c matrix,
  89 // or checks that a non-empty matrix is r×c. It zeroes
  90 // all the elements of the matrix.
  91 //
  92 // reuseAsZeroed must be kept in sync with reuseAs.
  93 func (m *Dense) reuseAsZeroed(r, c int) {
  94         if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
  95                 // Panic as a string, not a mat.Error.
  96                 panic("mat: caps not correctly set")
  97         }
  98         if m.IsZero() {
  99                 m.mat = blas64.General{
 100                         Rows:   r,
 101                         Cols:   c,
 102                         Stride: c,
 103                         Data:   useZeroed(m.mat.Data, r*c),
 104                 }
 105                 m.capRows = r
 106                 m.capCols = c
 107                 return
 108         }
 109         if r != m.mat.Rows || c != m.mat.Cols {
 110                 panic(ErrShape)
 111         }
 112         for i := 0; i < r; i++ {
 113                 zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c])
 114         }
 115 }
 116
 117 // untranspose untransposes a matrix if applicable. If a is an Untransposer, then
 118 // untranspose returns the underlying matrix and true. If it is not, then it returns
 119 // the input matrix and false.
 120 func untranspose(a Matrix) (Matrix, bool) {
 121         if ut, ok := a.(Untransposer); ok {
 122                 return ut.Untranspose(), true
 123         }
 124         return a, false
 125 }
 126
 127 // isolatedWorkspace returns a new dense matrix w with the size of a and
 128 // returns a callback to defer which performs cleanup at the return of the call.
 129 // This should be used when a method receiver is the same pointer as an input argument.
 130 func (m *Dense) isolatedWorkspace(a Matrix) (w *Dense, restore func()) {
 131         r, c := a.Dims()
 132         w = getWorkspace(r, c, false)
 133         return w, func() {
 134                 m.Copy(w)
 135                 putWorkspace(w)
 136         }
 137 }
 138
 139 // Reset zeros the dimensions of the matrix so that it can be reused as the
 140 // receiver of a dimensionally restricted operation.
 141 //
 142 // See the Reseter interface for more information.
 143 func (m *Dense) Reset() {
 144         // Row, Cols and Stride must be zeroed in unison.
 145         m.mat.Rows, m.mat.Cols, m.mat.Stride = 0, 0, 0
 146         m.capRows, m.capCols = 0, 0
 147         m.mat.Data = m.mat.Data[:0]
 148 }
 149
 150 // IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
 151 // receiver for size-restricted operations. Dense matrices can be zeroed using Reset.
 152 func (m *Dense) IsZero() bool {
 153         // It must be the case that m.Dims() returns
 154         // zeros in this case. See comment in Reset().
 155         return m.mat.Stride == 0
 156 }
 157
 158 // asTriDense returns a TriDense with the given size and side. The backing data
 159 // of the TriDense is the same as the receiver.
 160 func (m *Dense) asTriDense(n int, diag blas.Diag, uplo blas.Uplo) *TriDense {
 161         return &TriDense{
 162                 mat: blas64.Triangular{
 163                         N:      n,
 164                         Stride: m.mat.Stride,
 165                         Data:   m.mat.Data,
 166                         Uplo:   uplo,
 167                         Diag:   diag,
 168                 },
 169                 cap: n,
 170         }
 171 }
 172
 173 // DenseCopyOf returns a newly allocated copy of the elements of a.
 174 func DenseCopyOf(a Matrix) *Dense {
 175         d := &Dense{}
 176         d.Clone(a)
 177         return d
 178 }
 179
 180 // SetRawMatrix sets the underlying blas64.General used by the receiver.
 181 // Changes to elements in the receiver following the call will be reflected
 182 // in b.
 183 func (m *Dense) SetRawMatrix(b blas64.General) {
 184         m.capRows, m.capCols = b.Rows, b.Cols
 185         m.mat = b
 186 }
 187
 188 // RawMatrix returns the underlying blas64.General used by the receiver.
 189 // Changes to elements in the receiver following the call will be reflected
 190 // in returned blas64.General.
 191 func (m *Dense) RawMatrix() blas64.General { return m.mat }
 192
 193 // Dims returns the number of rows and columns in the matrix.
 194 func (m *Dense) Dims() (r, c int) { return m.mat.Rows, m.mat.Cols }
 195
 196 // Caps returns the number of rows and columns in the backing matrix.
 197 func (m *Dense) Caps() (r, c int) { return m.capRows, m.capCols }
 198
 199 // T performs an implicit transpose by returning the receiver inside a Transpose.
 200 func (m *Dense) T() Matrix {
 201         return Transpose{m}
 202 }
 203
 204 // ColView returns a Vector reflecting the column j, backed by the matrix data.
 205 //
 206 // See ColViewer for more information.
 207 func (m *Dense) ColView(j int) Vector {
 208         var v VecDense
 209         v.ColViewOf(m, j)
 210         return &v
 211 }
 212
 213 // SetCol sets the values in the specified column of the matrix to the values
 214 // in src. len(src) must equal the number of rows in the receiver.
 215 func (m *Dense) SetCol(j int, src []float64) {
 216         if j >= m.mat.Cols || j < 0 {
 217                 panic(ErrColAccess)
 218         }
 219         if len(src) != m.mat.Rows {
 220                 panic(ErrColLength)
 221         }
 222
 223         blas64.Copy(m.mat.Rows,
 224                 blas64.Vector{Inc: 1, Data: src},
 225                 blas64.Vector{Inc: m.mat.Stride, Data: m.mat.Data[j:]},
 226         )
 227 }
 228
 229 // SetRow sets the values in the specified rows of the matrix to the values
 230 // in src. len(src) must equal the number of columns in the receiver.
 231 func (m *Dense) SetRow(i int, src []float64) {
 232         if i >= m.mat.Rows || i < 0 {
 233                 panic(ErrRowAccess)
 234         }
 235         if len(src) != m.mat.Cols {
 236                 panic(ErrRowLength)
 237         }
 238
 239         copy(m.rawRowView(i), src)
 240 }
 241
 242 // RowView returns row i of the matrix data represented as a column vector,
 243 // backed by the matrix data.
 244 //
 245 // See RowViewer for more information.
 246 func (m *Dense) RowView(i int) Vector {
 247         var v VecDense
 248         v.RowViewOf(m, i)
 249         return &v
 250 }
 251
 252 // RawRowView returns a slice backed by the same array as backing the
 253 // receiver.
 254 func (m *Dense) RawRowView(i int) []float64 {
 255         if i >= m.mat.Rows || i < 0 {
 256                 panic(ErrRowAccess)
 257         }
 258         return m.rawRowView(i)
 259 }
 260
 261 func (m *Dense) rawRowView(i int) []float64 {
 262         return m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+m.mat.Cols]
 263 }
 264
 265 // Slice returns a new Matrix that shares backing data with the receiver.
 266 // The returned matrix starts at {i,j} of the receiver and extends k-i rows
 267 // and l-j columns. The final row in the resulting matrix is k-1 and the
 268 // final column is l-1.
 269 // Slice panics with ErrIndexOutOfRange if the slice is outside the capacity
 270 // of the receiver.
 271 func (m *Dense) Slice(i, k, j, l int) Matrix {
 272         mr, mc := m.Caps()
 273         if i < 0 || mr <= i || j < 0 || mc <= j || k <= i || mr < k || l <= j || mc < l {
 274                 panic(ErrIndexOutOfRange)
 275         }
 276         t := *m
 277         t.mat.Data = t.mat.Data[i*t.mat.Stride+j : (k-1)*t.mat.Stride+l]
 278         t.mat.Rows = k - i
 279         t.mat.Cols = l - j
 280         t.capRows -= i
 281         t.capCols -= j
 282         return &t
 283 }
 284
 285 // Grow returns the receiver expanded by r rows and c columns. If the dimensions
 286 // of the expanded matrix are outside the capacities of the receiver a new
 287 // allocation is made, otherwise not. Note the receiver itself is not modified
 288 // during the call to Grow.
 289 func (m *Dense) Grow(r, c int) Matrix {
 290         if r < 0 || c < 0 {
 291                 panic(ErrIndexOutOfRange)
 292         }
 293         if r == 0 && c == 0 {
 294                 return m
 295         }
 296
 297         r += m.mat.Rows
 298         c += m.mat.Cols
 299
 300         var t Dense
 301         switch {
 302         case m.mat.Rows == 0 || m.mat.Cols == 0:
 303                 t.mat = blas64.General{
 304                         Rows:   r,
 305                         Cols:   c,
 306                         Stride: c,
 307                         // We zero because we don't know how the matrix will be used.
 308                         // In other places, the mat is immediately filled with a result;
 309                         // this is not the case here.
 310                         Data: useZeroed(m.mat.Data, r*c),
 311                 }
 312         case r > m.capRows || c > m.capCols:
 313                 cr := max(r, m.capRows)
 314                 cc := max(c, m.capCols)
 315                 t.mat = blas64.General{
 316                         Rows:   r,
 317                         Cols:   c,
 318                         Stride: cc,
 319                         Data:   make([]float64, cr*cc),
 320                 }
 321                 t.capRows = cr
 322                 t.capCols = cc
 323                 // Copy the complete matrix over to the new matrix.
 324                 // Including elements not currently visible. Use a temporary structure
 325                 // to avoid modifying the receiver.
 326                 var tmp Dense
 327                 tmp.mat = blas64.General{
 328                         Rows:   m.mat.Rows,
 329                         Cols:   m.mat.Cols,
 330                         Stride: m.mat.Stride,
 331                         Data:   m.mat.Data,
 332                 }
 333                 tmp.capRows = m.capRows
 334                 tmp.capCols = m.capCols
 335                 t.Copy(&tmp)
 336                 return &t
 337         default:
 338                 t.mat = blas64.General{
 339                         Data:   m.mat.Data[:(r-1)*m.mat.Stride+c],
 340                         Rows:   r,
 341                         Cols:   c,
 342                         Stride: m.mat.Stride,
 343                 }
 344         }
 345         t.capRows = r
 346         t.capCols = c
 347         return &t
 348 }
 349
 350 // Clone makes a copy of a into the receiver, overwriting the previous value of
 351 // the receiver. The clone operation does not make any restriction on shape and
 352 // will not cause shadowing.
 353 //
 354 // See the Cloner interface for more information.
 355 func (m *Dense) Clone(a Matrix) {
 356         r, c := a.Dims()
 357         mat := blas64.General{
 358                 Rows:   r,
 359                 Cols:   c,
 360                 Stride: c,
 361         }
 362         m.capRows, m.capCols = r, c
 363
 364         aU, trans := untranspose(a)
 365         switch aU := aU.(type) {
 366         case RawMatrixer:
 367                 amat := aU.RawMatrix()
 368                 mat.Data = make([]float64, r*c)
 369                 if trans {
 370                         for i := 0; i < r; i++ {
 371                                 blas64.Copy(c,
 372                                         blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
 373                                         blas64.Vector{Inc: 1, Data: mat.Data[i*c : (i+1)*c]})
 374                         }
 375                 } else {
 376                         for i := 0; i < r; i++ {
 377                                 copy(mat.Data[i*c:(i+1)*c], amat.Data[i*amat.Stride:i*amat.Stride+c])
 378                         }
 379                 }
 380         case *VecDense:
 381                 amat := aU.mat
 382                 mat.Data = make([]float64, aU.n)
 383                 blas64.Copy(aU.n,
 384                         blas64.Vector{Inc: amat.Inc, Data: amat.Data},
 385                         blas64.Vector{Inc: 1, Data: mat.Data})
 386         default:
 387                 mat.Data = make([]float64, r*c)
 388                 w := *m
 389                 w.mat = mat
 390                 for i := 0; i < r; i++ {
 391                         for j := 0; j < c; j++ {
 392                                 w.set(i, j, a.At(i, j))
 393                         }
 394                 }
 395                 *m = w
 396                 return
 397         }
 398         m.mat = mat
 399 }
 400
 401 // Copy makes a copy of elements of a into the receiver. It is similar to the
 402 // built-in copy; it copies as much as the overlap between the two matrices and
 403 // returns the number of rows and columns it copied. If a aliases the receiver
 404 // and is a transposed Dense or VecDense, with a non-unitary increment, Copy will
 405 // panic.
 406 //
 407 // See the Copier interface for more information.
 408 func (m *Dense) Copy(a Matrix) (r, c int) {
 409         r, c = a.Dims()
 410         if a == m {
 411                 return r, c
 412         }
 413         r = min(r, m.mat.Rows)
 414         c = min(c, m.mat.Cols)
 415         if r == 0 || c == 0 {
 416                 return 0, 0
 417         }
 418
 419         aU, trans := untranspose(a)
 420         switch aU := aU.(type) {
 421         case RawMatrixer:
 422                 amat := aU.RawMatrix()
 423                 if trans {
 424                         if amat.Stride != 1 {
 425                                 m.checkOverlap(amat)
 426                         }
 427                         for i := 0; i < r; i++ {
 428                                 blas64.Copy(c,
 429                                         blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
 430                                         blas64.Vector{Inc: 1, Data: m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c]})
 431                         }
 432                 } else {
 433                         switch o := offset(m.mat.Data, amat.Data); {
 434                         case o < 0:
 435                                 for i := r - 1; i >= 0; i-- {
 436                                         copy(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], amat.Data[i*amat.Stride:i*amat.Stride+c])
 437                                 }
 438                         case o > 0:
 439                                 for i := 0; i < r; i++ {
 440                                         copy(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], amat.Data[i*amat.Stride:i*amat.Stride+c])
 441                                 }
 442                         default:
 443                                 // Nothing to do.
 444                         }
 445                 }
 446         case *VecDense:
 447                 var n, stride int
 448                 amat := aU.mat
 449                 if trans {
 450                         if amat.Inc != 1 {
 451                                 m.checkOverlap(aU.asGeneral())
 452                         }
 453                         n = c
 454                         stride = 1
 455                 } else {
 456                         n = r
 457                         stride = m.mat.Stride
 458                 }
 459                 if amat.Inc == 1 && stride == 1 {
 460                         copy(m.mat.Data, amat.Data[:n])
 461                         break
 462                 }
 463                 switch o := offset(m.mat.Data, amat.Data); {
 464                 case o < 0:
 465                         blas64.Copy(n,
 466                                 blas64.Vector{Inc: -amat.Inc, Data: amat.Data},
 467                                 blas64.Vector{Inc: -stride, Data: m.mat.Data})
 468                 case o > 0:
 469                         blas64.Copy(n,
 470                                 blas64.Vector{Inc: amat.Inc, Data: amat.Data},
 471                                 blas64.Vector{Inc: stride, Data: m.mat.Data})
 472                 default:
 473                         // Nothing to do.
 474                 }
 475         default:
 476                 for i := 0; i < r; i++ {
 477                         for j := 0; j < c; j++ {
 478                                 m.set(i, j, a.At(i, j))
 479                         }
 480                 }
 481         }
 482
 483         return r, c
 484 }
 485
 486 // Stack appends the rows of b onto the rows of a, placing the result into the
 487 // receiver with b placed in the greater indexed rows. Stack will panic if the
 488 // two input matrices do not have the same number of columns or the constructed
 489 // stacked matrix is not the same shape as the receiver.
 490 func (m *Dense) Stack(a, b Matrix) {
 491         ar, ac := a.Dims()
 492         br, bc := b.Dims()
 493         if ac != bc || m == a || m == b {
 494                 panic(ErrShape)
 495         }
 496
 497         m.reuseAs(ar+br, ac)
 498
 499         m.Copy(a)
 500         w := m.Slice(ar, ar+br, 0, bc).(*Dense)
 501         w.Copy(b)
 502 }
 503
 504 // Augment creates the augmented matrix of a and b, where b is placed in the
 505 // greater indexed columns. Augment will panic if the two input matrices do
 506 // not have the same number of rows or the constructed augmented matrix is
 507 // not the same shape as the receiver.
 508 func (m *Dense) Augment(a, b Matrix) {
 509         ar, ac := a.Dims()
 510         br, bc := b.Dims()
 511         if ar != br || m == a || m == b {
 512                 panic(ErrShape)
 513         }
 514
 515         m.reuseAs(ar, ac+bc)
 516
 517         m.Copy(a)
 518         w := m.Slice(0, br, ac, ac+bc).(*Dense)
 519         w.Copy(b)
 520 }