triband.go 11 KB

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