band.go 7.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264
  1. // Copyright ©2017 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/blas64"
  7. )
  8. var (
  9. bandDense *BandDense
  10. _ Matrix = bandDense
  11. _ Banded = bandDense
  12. _ RawBander = bandDense
  13. _ NonZeroDoer = bandDense
  14. _ RowNonZeroDoer = bandDense
  15. _ ColNonZeroDoer = bandDense
  16. )
  17. // BandDense represents a band matrix in dense storage format.
  18. type BandDense struct {
  19. mat blas64.Band
  20. }
  21. // Banded is a band matrix representation.
  22. type Banded interface {
  23. Matrix
  24. // Bandwidth returns the lower and upper bandwidth values for
  25. // the matrix. The total bandwidth of the matrix is kl+ku+1.
  26. Bandwidth() (kl, ku int)
  27. // TBand is the equivalent of the T() method in the Matrix
  28. // interface but guarantees the transpose is of banded type.
  29. TBand() Banded
  30. }
  31. // A RawBander can return a blas64.Band representation of the receiver.
  32. // Changes to the blas64.Band.Data slice will be reflected in the original
  33. // matrix, changes to the Rows, Cols, KL, KU and Stride fields will not.
  34. type RawBander interface {
  35. RawBand() blas64.Band
  36. }
  37. // A MutableBanded can set elements of a band matrix.
  38. type MutableBanded interface {
  39. Banded
  40. SetBand(i, j int, v float64)
  41. }
  42. var (
  43. _ Matrix = TransposeBand{}
  44. _ Banded = TransposeBand{}
  45. _ UntransposeBander = TransposeBand{}
  46. )
  47. // TransposeBand is a type for performing an implicit transpose of a band
  48. // matrix. It implements the Banded interface, returning values from the
  49. // transpose of the matrix within.
  50. type TransposeBand struct {
  51. Banded Banded
  52. }
  53. // At returns the value of the element at row i and column j of the transposed
  54. // matrix, that is, row j and column i of the Banded field.
  55. func (t TransposeBand) At(i, j int) float64 {
  56. return t.Banded.At(j, i)
  57. }
  58. // Dims returns the dimensions of the transposed matrix.
  59. func (t TransposeBand) Dims() (r, c int) {
  60. c, r = t.Banded.Dims()
  61. return r, c
  62. }
  63. // T performs an implicit transpose by returning the Banded field.
  64. func (t TransposeBand) T() Matrix {
  65. return t.Banded
  66. }
  67. // Bandwidth returns the lower and upper bandwidth values for
  68. // the transposed matrix.
  69. func (t TransposeBand) Bandwidth() (kl, ku int) {
  70. kl, ku = t.Banded.Bandwidth()
  71. return ku, kl
  72. }
  73. // TBand performs an implicit transpose by returning the Banded field.
  74. func (t TransposeBand) TBand() Banded {
  75. return t.Banded
  76. }
  77. // Untranspose returns the Banded field.
  78. func (t TransposeBand) Untranspose() Matrix {
  79. return t.Banded
  80. }
  81. // UntransposeBand returns the Banded field.
  82. func (t TransposeBand) UntransposeBand() Banded {
  83. return t.Banded
  84. }
  85. // NewBandDense creates a new Band matrix with r rows and c columns. If data == nil,
  86. // a new slice is allocated for the backing slice. If len(data) == min(r, c+kl)*(kl+ku+1),
  87. // data is used as the backing slice, and changes to the elements of the returned
  88. // BandDense will be reflected in data. If neither of these is true, NewBandDense
  89. // will panic. kl must be at least zero and less r, and ku must be at least zero and
  90. // less than c, otherwise NewBandDense will panic.
  91. // NewBandDense will panic if either r or c is zero.
  92. //
  93. // The data must be arranged in row-major order constructed by removing the zeros
  94. // from the rows outside the band and aligning the diagonals. For example, the matrix
  95. // 1 2 3 0 0 0
  96. // 4 5 6 7 0 0
  97. // 0 8 9 10 11 0
  98. // 0 0 12 13 14 15
  99. // 0 0 0 16 17 18
  100. // 0 0 0 0 19 20
  101. // becomes (* entries are never accessed)
  102. // * 1 2 3
  103. // 4 5 6 7
  104. // 8 9 10 11
  105. // 12 13 14 15
  106. // 16 17 18 *
  107. // 19 20 * *
  108. // which is passed to NewBandDense as []float64{*, 1, 2, 3, 4, ...} with kl=1 and ku=2.
  109. // Only the values in the band portion of the matrix are used.
  110. func NewBandDense(r, c, kl, ku int, data []float64) *BandDense {
  111. if r <= 0 || c <= 0 || kl < 0 || ku < 0 {
  112. if r == 0 || c == 0 {
  113. panic(ErrZeroLength)
  114. }
  115. panic("mat: negative dimension")
  116. }
  117. if kl+1 > r || ku+1 > c {
  118. panic("mat: band out of range")
  119. }
  120. bc := kl + ku + 1
  121. if data != nil && len(data) != min(r, c+kl)*bc {
  122. panic(ErrShape)
  123. }
  124. if data == nil {
  125. data = make([]float64, min(r, c+kl)*bc)
  126. }
  127. return &BandDense{
  128. mat: blas64.Band{
  129. Rows: r,
  130. Cols: c,
  131. KL: kl,
  132. KU: ku,
  133. Stride: bc,
  134. Data: data,
  135. },
  136. }
  137. }
  138. // NewDiagonalRect is a convenience function that returns a diagonal matrix represented by a
  139. // BandDense. The length of data must be min(r, c) otherwise NewDiagonalRect will panic.
  140. func NewDiagonalRect(r, c int, data []float64) *BandDense {
  141. return NewBandDense(r, c, 0, 0, data)
  142. }
  143. // Dims returns the number of rows and columns in the matrix.
  144. func (b *BandDense) Dims() (r, c int) {
  145. return b.mat.Rows, b.mat.Cols
  146. }
  147. // Bandwidth returns the upper and lower bandwidths of the matrix.
  148. func (b *BandDense) Bandwidth() (kl, ku int) {
  149. return b.mat.KL, b.mat.KU
  150. }
  151. // T performs an implicit transpose by returning the receiver inside a Transpose.
  152. func (b *BandDense) T() Matrix {
  153. return Transpose{b}
  154. }
  155. // TBand performs an implicit transpose by returning the receiver inside a TransposeBand.
  156. func (b *BandDense) TBand() Banded {
  157. return TransposeBand{b}
  158. }
  159. // RawBand returns the underlying blas64.Band used by the receiver.
  160. // Changes to elements in the receiver following the call will be reflected
  161. // in returned blas64.Band.
  162. func (b *BandDense) RawBand() blas64.Band {
  163. return b.mat
  164. }
  165. // SetRawBand sets the underlying blas64.Band used by the receiver.
  166. // Changes to elements in the receiver following the call will be reflected
  167. // in the input.
  168. func (b *BandDense) SetRawBand(mat blas64.Band) {
  169. b.mat = mat
  170. }
  171. // DiagView returns the diagonal as a matrix backed by the original data.
  172. func (b *BandDense) DiagView() Diagonal {
  173. n := min(b.mat.Rows, b.mat.Cols)
  174. return &DiagDense{
  175. mat: blas64.Vector{
  176. N: n,
  177. Inc: b.mat.Stride,
  178. Data: b.mat.Data[b.mat.KL : (n-1)*b.mat.Stride+b.mat.KL+1],
  179. },
  180. }
  181. }
  182. // DoNonZero calls the function fn for each of the non-zero elements of b. The function fn
  183. // takes a row/column index and the element value of b at (i, j).
  184. func (b *BandDense) DoNonZero(fn func(i, j int, v float64)) {
  185. for i := 0; i < min(b.mat.Rows, b.mat.Cols+b.mat.KL); i++ {
  186. for j := max(0, i-b.mat.KL); j < min(b.mat.Cols, i+b.mat.KU+1); j++ {
  187. v := b.at(i, j)
  188. if v != 0 {
  189. fn(i, j, v)
  190. }
  191. }
  192. }
  193. }
  194. // DoRowNonZero calls the function fn for each of the non-zero elements of row i of b. The function fn
  195. // takes a row/column index and the element value of b at (i, j).
  196. func (b *BandDense) DoRowNonZero(i int, fn func(i, j int, v float64)) {
  197. if i < 0 || b.mat.Rows <= i {
  198. panic(ErrRowAccess)
  199. }
  200. for j := max(0, i-b.mat.KL); j < min(b.mat.Cols, i+b.mat.KU+1); j++ {
  201. v := b.at(i, j)
  202. if v != 0 {
  203. fn(i, j, v)
  204. }
  205. }
  206. }
  207. // DoColNonZero calls the function fn for each of the non-zero elements of column j of b. The function fn
  208. // takes a row/column index and the element value of b at (i, j).
  209. func (b *BandDense) DoColNonZero(j int, fn func(i, j int, v float64)) {
  210. if j < 0 || b.mat.Cols <= j {
  211. panic(ErrColAccess)
  212. }
  213. for i := 0; i < min(b.mat.Rows, b.mat.Cols+b.mat.KL); i++ {
  214. if i-b.mat.KL <= j && j < i+b.mat.KU+1 {
  215. v := b.at(i, j)
  216. if v != 0 {
  217. fn(i, j, v)
  218. }
  219. }
  220. }
  221. }
  222. // Zero sets all of the matrix elements to zero.
  223. func (b *BandDense) Zero() {
  224. m := b.mat.Rows
  225. kL := b.mat.KL
  226. nCol := b.mat.KU + 1 + kL
  227. for i := 0; i < m; i++ {
  228. l := max(0, kL-i)
  229. u := min(nCol, m+kL-i)
  230. zero(b.mat.Data[i*b.mat.Stride+l : i*b.mat.Stride+u])
  231. }
  232. }