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