dense.go 15 KB

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