dense.go 14 KB

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