triband.go 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354
  1. // Copyright ©2018 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. triBand TriBanded
  11. _ Banded = triBand
  12. _ Triangular = triBand
  13. triBandDense *TriBandDense
  14. _ Matrix = triBandDense
  15. _ Triangular = triBandDense
  16. _ Banded = triBandDense
  17. _ TriBanded = triBandDense
  18. _ RawTriBander = triBandDense
  19. _ MutableTriBanded = triBandDense
  20. )
  21. // TriBanded is a triangular band matrix interface type.
  22. type TriBanded interface {
  23. Banded
  24. // Triangle returns the number of rows/columns in the matrix and its
  25. // orientation.
  26. Triangle() (n int, kind TriKind)
  27. // TTri is the equivalent of the T() method in the Matrix interface but
  28. // guarantees the transpose is of triangular type.
  29. TTri() Triangular
  30. // TriBand returns the number of rows/columns in the matrix, the
  31. // size of the bandwidth, and the orientation.
  32. TriBand() (n, k int, kind TriKind)
  33. // TTriBand is the equivalent of the T() method in the Matrix interface but
  34. // guarantees the transpose is of banded triangular type.
  35. TTriBand() TriBanded
  36. }
  37. // A RawTriBander can return a blas64.TriangularBand representation of the receiver.
  38. // Changes to the blas64.TriangularBand.Data slice will be reflected in the original
  39. // matrix, changes to the N, K, Stride, Uplo and Diag fields will not.
  40. type RawTriBander interface {
  41. RawTriBand() blas64.TriangularBand
  42. }
  43. // MutableTriBanded is a triangular band matrix interface type that allows
  44. // elements to be altered.
  45. type MutableTriBanded interface {
  46. TriBanded
  47. SetTriBand(i, j int, v float64)
  48. }
  49. var (
  50. tTriBand TransposeTriBand
  51. _ Matrix = tTriBand
  52. _ TriBanded = tTriBand
  53. _ Untransposer = tTriBand
  54. _ UntransposeTrier = tTriBand
  55. _ UntransposeBander = tTriBand
  56. _ UntransposeTriBander = tTriBand
  57. )
  58. // TransposeTriBand is a type for performing an implicit transpose of a TriBanded
  59. // matrix. It implements the TriBanded interface, returning values from the
  60. // transpose of the matrix within.
  61. type TransposeTriBand struct {
  62. TriBanded TriBanded
  63. }
  64. // At returns the value of the element at row i and column j of the transposed
  65. // matrix, that is, row j and column i of the TriBanded field.
  66. func (t TransposeTriBand) At(i, j int) float64 {
  67. return t.TriBanded.At(j, i)
  68. }
  69. // Dims returns the dimensions of the transposed matrix. TriBanded matrices are
  70. // square and thus this is the same size as the original TriBanded.
  71. func (t TransposeTriBand) Dims() (r, c int) {
  72. c, r = t.TriBanded.Dims()
  73. return r, c
  74. }
  75. // T performs an implicit transpose by returning the TriBand field.
  76. func (t TransposeTriBand) T() Matrix {
  77. return t.TriBanded
  78. }
  79. // Triangle returns the number of rows/columns in the matrix and its orientation.
  80. func (t TransposeTriBand) Triangle() (int, TriKind) {
  81. n, upper := t.TriBanded.Triangle()
  82. return n, !upper
  83. }
  84. // TTri performs an implicit transpose by returning the TriBand field.
  85. func (t TransposeTriBand) TTri() Triangular {
  86. return t.TriBanded
  87. }
  88. // Bandwidth returns the upper and lower bandwidths of the matrix.
  89. func (t TransposeTriBand) Bandwidth() (kl, ku int) {
  90. kl, ku = t.TriBanded.Bandwidth()
  91. return ku, kl
  92. }
  93. // TBand performs an implicit transpose by returning the TriBand field.
  94. func (t TransposeTriBand) TBand() Banded {
  95. return t.TriBanded
  96. }
  97. // TriBand returns the number of rows/columns in the matrix, the
  98. // size of the bandwidth, and the orientation.
  99. func (t TransposeTriBand) TriBand() (n, k int, kind TriKind) {
  100. n, k, kind = t.TriBanded.TriBand()
  101. return n, k, !kind
  102. }
  103. // TTriBand performs an implicit transpose by returning the TriBand field.
  104. func (t TransposeTriBand) TTriBand() TriBanded {
  105. return t.TriBanded
  106. }
  107. // Untranspose returns the Triangular field.
  108. func (t TransposeTriBand) Untranspose() Matrix {
  109. return t.TriBanded
  110. }
  111. // UntransposeTri returns the underlying Triangular matrix.
  112. func (t TransposeTriBand) UntransposeTri() Triangular {
  113. return t.TriBanded
  114. }
  115. // UntransposeBand returns the underlying Banded matrix.
  116. func (t TransposeTriBand) UntransposeBand() Banded {
  117. return t.TriBanded
  118. }
  119. // UntransposeTriBand returns the underlying TriBanded matrix.
  120. func (t TransposeTriBand) UntransposeTriBand() TriBanded {
  121. return t.TriBanded
  122. }
  123. // TriBandDense represents a triangular band matrix in dense storage format.
  124. type TriBandDense struct {
  125. mat blas64.TriangularBand
  126. }
  127. // NewTriBandDense creates a new triangular banded matrix with n rows and columns,
  128. // k bands in the direction of the specified kind. If data == nil,
  129. // a new slice is allocated for the backing slice. If len(data) == n*(k+1),
  130. // data is used as the backing slice, and changes to the elements of the returned
  131. // TriBandDense will be reflected in data. If neither of these is true, NewTriBandDense
  132. // will panic. k must be at least zero and less than n, otherwise NewTriBandDense will panic.
  133. //
  134. // The data must be arranged in row-major order constructed by removing the zeros
  135. // from the rows outside the band and aligning the diagonals. For example, if
  136. // the upper-triangular banded matrix
  137. // 1 2 3 0 0 0
  138. // 0 4 5 6 0 0
  139. // 0 0 7 8 9 0
  140. // 0 0 0 10 11 12
  141. // 0 0 0 0 13 14
  142. // 0 0 0 0 0 15
  143. // becomes (* entries are never accessed)
  144. // 1 2 3
  145. // 4 5 6
  146. // 7 8 9
  147. // 10 11 12
  148. // 13 14 *
  149. // 15 * *
  150. // which is passed to NewTriBandDense as []float64{1, 2, ..., 15, *, *, *}
  151. // with k=2 and kind = mat.Upper.
  152. // The lower triangular banded matrix
  153. // 1 0 0 0 0 0
  154. // 2 3 0 0 0 0
  155. // 4 5 6 0 0 0
  156. // 0 7 8 9 0 0
  157. // 0 0 10 11 12 0
  158. // 0 0 0 13 14 15
  159. // becomes (* entries are never accessed)
  160. // * * 1
  161. // * 2 3
  162. // 4 5 6
  163. // 7 8 9
  164. // 10 11 12
  165. // 13 14 15
  166. // which is passed to NewTriBandDense as []float64{*, *, *, 1, 2, ..., 15}
  167. // with k=2 and kind = mat.Lower.
  168. // Only the values in the band portion of the matrix are used.
  169. func NewTriBandDense(n, k int, kind TriKind, data []float64) *TriBandDense {
  170. if n <= 0 || k < 0 {
  171. if n == 0 {
  172. panic(ErrZeroLength)
  173. }
  174. panic("mat: negative dimension")
  175. }
  176. if k+1 > n {
  177. panic("mat: band out of range")
  178. }
  179. bc := k + 1
  180. if data != nil && len(data) != n*bc {
  181. panic(ErrShape)
  182. }
  183. if data == nil {
  184. data = make([]float64, n*bc)
  185. }
  186. uplo := blas.Lower
  187. if kind {
  188. uplo = blas.Upper
  189. }
  190. return &TriBandDense{
  191. mat: blas64.TriangularBand{
  192. Uplo: uplo,
  193. Diag: blas.NonUnit,
  194. N: n,
  195. K: k,
  196. Data: data,
  197. Stride: bc,
  198. },
  199. }
  200. }
  201. // Dims returns the number of rows and columns in the matrix.
  202. func (t *TriBandDense) Dims() (r, c int) {
  203. return t.mat.N, t.mat.N
  204. }
  205. // T performs an implicit transpose by returning the receiver inside a Transpose.
  206. func (t *TriBandDense) T() Matrix {
  207. return Transpose{t}
  208. }
  209. // IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
  210. // receiver for size-restricted operations. TriBandDense matrices can be zeroed using Reset.
  211. func (t *TriBandDense) IsZero() bool {
  212. // It must be the case that t.Dims() returns
  213. // zeros in this case. See comment in Reset().
  214. return t.mat.Stride == 0
  215. }
  216. // Reset zeros the dimensions of the matrix so that it can be reused as the
  217. // receiver of a dimensionally restricted operation.
  218. //
  219. // See the Reseter interface for more information.
  220. func (t *TriBandDense) Reset() {
  221. t.mat.N = 0
  222. t.mat.Stride = 0
  223. t.mat.K = 0
  224. t.mat.Data = t.mat.Data[:0]
  225. }
  226. // Zero sets all of the matrix elements to zero.
  227. func (t *TriBandDense) Zero() {
  228. if t.isUpper() {
  229. for i := 0; i < t.mat.N; i++ {
  230. u := min(1+t.mat.K, t.mat.N-i)
  231. zero(t.mat.Data[i*t.mat.Stride : i*t.mat.Stride+u])
  232. }
  233. return
  234. }
  235. for i := 0; i < t.mat.N; i++ {
  236. l := max(0, t.mat.K-i)
  237. zero(t.mat.Data[i*t.mat.Stride+l : i*t.mat.Stride+t.mat.K+1])
  238. }
  239. }
  240. func (t *TriBandDense) isUpper() bool {
  241. return isUpperUplo(t.mat.Uplo)
  242. }
  243. func (t *TriBandDense) triKind() TriKind {
  244. return TriKind(isUpperUplo(t.mat.Uplo))
  245. }
  246. // Triangle returns the dimension of t and its orientation. The returned
  247. // orientation is only valid when n is not zero.
  248. func (t *TriBandDense) Triangle() (n int, kind TriKind) {
  249. return t.mat.N, t.triKind()
  250. }
  251. // TTri performs an implicit transpose by returning the receiver inside a TransposeTri.
  252. func (t *TriBandDense) TTri() Triangular {
  253. return TransposeTri{t}
  254. }
  255. // Bandwidth returns the upper and lower bandwidths of the matrix.
  256. func (t *TriBandDense) Bandwidth() (kl, ku int) {
  257. if t.isUpper() {
  258. return 0, t.mat.K
  259. }
  260. return t.mat.K, 0
  261. }
  262. // TBand performs an implicit transpose by returning the receiver inside a TransposeBand.
  263. func (t *TriBandDense) TBand() Banded {
  264. return TransposeBand{t}
  265. }
  266. // TriBand returns the number of rows/columns in the matrix, the
  267. // size of the bandwidth, and the orientation.
  268. func (t *TriBandDense) TriBand() (n, k int, kind TriKind) {
  269. return t.mat.N, t.mat.K, TriKind(!t.IsZero()) && t.triKind()
  270. }
  271. // TTriBand performs an implicit transpose by returning the receiver inside a TransposeTriBand.
  272. func (t *TriBandDense) TTriBand() TriBanded {
  273. return TransposeTriBand{t}
  274. }
  275. // RawTriBand returns the underlying blas64.TriangularBand used by the receiver.
  276. // Changes to the blas64.TriangularBand.Data slice will be reflected in the original
  277. // matrix, changes to the N, K, Stride, Uplo and Diag fields will not.
  278. func (t *TriBandDense) RawTriBand() blas64.TriangularBand {
  279. return t.mat
  280. }
  281. // SetRawTriBand sets the underlying blas64.TriangularBand used by the receiver.
  282. // Changes to elements in the receiver following the call will be reflected
  283. // in the input.
  284. //
  285. // The supplied TriangularBand must not use blas.Unit storage format.
  286. func (t *TriBandDense) SetRawTriBand(mat blas64.TriangularBand) {
  287. if mat.Diag == blas.Unit {
  288. panic("mat: cannot set TriBand with Unit storage")
  289. }
  290. t.mat = mat
  291. }
  292. // DiagView returns the diagonal as a matrix backed by the original data.
  293. func (t *TriBandDense) DiagView() Diagonal {
  294. if t.mat.Diag == blas.Unit {
  295. panic("mat: cannot take view of Unit diagonal")
  296. }
  297. n := t.mat.N
  298. data := t.mat.Data
  299. if !t.isUpper() {
  300. data = data[t.mat.K:]
  301. }
  302. return &DiagDense{
  303. mat: blas64.Vector{
  304. N: n,
  305. Inc: t.mat.Stride,
  306. Data: data[:(n-1)*t.mat.Stride+1],
  307. },
  308. }
  309. }