vendor/gonum.org/v1/gonum/mat/triangular.go

   1 // Copyright ©2015 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
   5 package mat
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/blas/blas64"
  12         "gonum.org/v1/gonum/lapack/lapack64"
  13 )
  14
  15 var (
  16         triDense *TriDense
  17         _        Matrix            = triDense
  18         _        Triangular        = triDense
  19         _        RawTriangular     = triDense
  20         _        MutableTriangular = triDense
  21
  22         _ NonZeroDoer    = triDense
  23         _ RowNonZeroDoer = triDense
  24         _ ColNonZeroDoer = triDense
  25 )
  26
  27 const badTriCap = "mat: bad capacity for TriDense"
  28
  29 // TriDense represents an upper or lower triangular matrix in dense storage
  30 // format.
  31 type TriDense struct {
  32         mat blas64.Triangular
  33         cap int
  34 }
  35
  36 // Triangular represents a triangular matrix. Triangular matrices are always square.
  37 type Triangular interface {
  38         Matrix
  39         // Triangular returns the number of rows/columns in the matrix and its
  40         // orientation.
  41         Triangle() (n int, kind TriKind)
  42
  43         // TTri is the equivalent of the T() method in the Matrix interface but
  44         // guarantees the transpose is of triangular type.
  45         TTri() Triangular
  46 }
  47
  48 // A RawTriangular can return a view of itself as a BLAS Triangular matrix.
  49 type RawTriangular interface {
  50         RawTriangular() blas64.Triangular
  51 }
  52
  53 // A MutableTriangular can set elements of a triangular matrix.
  54 type MutableTriangular interface {
  55         Triangular
  56         SetTri(i, j int, v float64)
  57 }
  58
  59 var (
  60         _ Matrix           = TransposeTri{}
  61         _ Triangular       = TransposeTri{}
  62         _ UntransposeTrier = TransposeTri{}
  63 )
  64
  65 // TransposeTri is a type for performing an implicit transpose of a Triangular
  66 // matrix. It implements the Triangular interface, returning values from the
  67 // transpose of the matrix within.
  68 type TransposeTri struct {
  69         Triangular Triangular
  70 }
  71
  72 // At returns the value of the element at row i and column j of the transposed
  73 // matrix, that is, row j and column i of the Triangular field.
  74 func (t TransposeTri) At(i, j int) float64 {
  75         return t.Triangular.At(j, i)
  76 }
  77
  78 // Dims returns the dimensions of the transposed matrix. Triangular matrices are
  79 // square and thus this is the same size as the original Triangular.
  80 func (t TransposeTri) Dims() (r, c int) {
  81         c, r = t.Triangular.Dims()
  82         return r, c
  83 }
  84
  85 // T performs an implicit transpose by returning the Triangular field.
  86 func (t TransposeTri) T() Matrix {
  87         return t.Triangular
  88 }
  89
  90 // Triangle returns the number of rows/columns in the matrix and its orientation.
  91 func (t TransposeTri) Triangle() (int, TriKind) {
  92         n, upper := t.Triangular.Triangle()
  93         return n, !upper
  94 }
  95
  96 // TTri performs an implicit transpose by returning the Triangular field.
  97 func (t TransposeTri) TTri() Triangular {
  98         return t.Triangular
  99 }
 100
 101 // Untranspose returns the Triangular field.
 102 func (t TransposeTri) Untranspose() Matrix {
 103         return t.Triangular
 104 }
 105
 106 func (t TransposeTri) UntransposeTri() Triangular {
 107         return t.Triangular
 108 }
 109
 110 // NewTriDense creates a new Triangular matrix with n rows and columns. If data == nil,
 111 // a new slice is allocated for the backing slice. If len(data) == n*n, data is
 112 // used as the backing slice, and changes to the elements of the returned TriDense
 113 // will be reflected in data. If neither of these is true, NewTriDense will panic.
 114 //
 115 // The data must be arranged in row-major order, i.e. the (i*c + j)-th
 116 // element in the data slice is the {i, j}-th element in the matrix.
 117 // Only the values in the triangular portion corresponding to kind are used.
 118 func NewTriDense(n int, kind TriKind, data []float64) *TriDense {
 119         if n < 0 {
 120                 panic("mat: negative dimension")
 121         }
 122         if data != nil && len(data) != n*n {
 123                 panic(ErrShape)
 124         }
 125         if data == nil {
 126                 data = make([]float64, n*n)
 127         }
 128         uplo := blas.Lower
 129         if kind == Upper {
 130                 uplo = blas.Upper
 131         }
 132         return &TriDense{
 133                 mat: blas64.Triangular{
 134                         N:      n,
 135                         Stride: n,
 136                         Data:   data,
 137                         Uplo:   uplo,
 138                         Diag:   blas.NonUnit,
 139                 },
 140                 cap: n,
 141         }
 142 }
 143
 144 func (t *TriDense) Dims() (r, c int) {
 145         return t.mat.N, t.mat.N
 146 }
 147
 148 // Triangle returns the dimension of t and its orientation. The returned
 149 // orientation is only valid when n is not zero.
 150 func (t *TriDense) Triangle() (n int, kind TriKind) {
 151         return t.mat.N, TriKind(!t.IsZero()) && t.triKind()
 152 }
 153
 154 func (t *TriDense) isUpper() bool {
 155         return isUpperUplo(t.mat.Uplo)
 156 }
 157
 158 func (t *TriDense) triKind() TriKind {
 159         return TriKind(isUpperUplo(t.mat.Uplo))
 160 }
 161
 162 func isUpperUplo(u blas.Uplo) bool {
 163         switch u {
 164         case blas.Upper:
 165                 return true
 166         case blas.Lower:
 167                 return false
 168         default:
 169                 panic(badTriangle)
 170         }
 171 }
 172
 173 // asSymBlas returns the receiver restructured as a blas64.Symmetric with the
 174 // same backing memory. Panics if the receiver is unit.
 175 // This returns a blas64.Symmetric and not a *SymDense because SymDense can only
 176 // be upper triangular.
 177 func (t *TriDense) asSymBlas() blas64.Symmetric {
 178         if t.mat.Diag == blas.Unit {
 179                 panic("mat: cannot convert unit TriDense into blas64.Symmetric")
 180         }
 181         return blas64.Symmetric{
 182                 N:      t.mat.N,
 183                 Stride: t.mat.Stride,
 184                 Data:   t.mat.Data,
 185                 Uplo:   t.mat.Uplo,
 186         }
 187 }
 188
 189 // T performs an implicit transpose by returning the receiver inside a Transpose.
 190 func (t *TriDense) T() Matrix {
 191         return Transpose{t}
 192 }
 193
 194 // TTri performs an implicit transpose by returning the receiver inside a TransposeTri.
 195 func (t *TriDense) TTri() Triangular {
 196         return TransposeTri{t}
 197 }
 198
 199 func (t *TriDense) RawTriangular() blas64.Triangular {
 200         return t.mat
 201 }
 202
 203 // Reset zeros the dimensions of the matrix so that it can be reused as the
 204 // receiver of a dimensionally restricted operation.
 205 //
 206 // See the Reseter interface for more information.
 207 func (t *TriDense) Reset() {
 208         // N and Stride must be zeroed in unison.
 209         t.mat.N, t.mat.Stride = 0, 0
 210         // Defensively zero Uplo to ensure
 211         // it is set correctly later.
 212         t.mat.Uplo = 0
 213         t.mat.Data = t.mat.Data[:0]
 214 }
 215
 216 // IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
 217 // receiver for size-restricted operations. TriDense matrices can be zeroed using Reset.
 218 func (t *TriDense) IsZero() bool {
 219         // It must be the case that t.Dims() returns
 220         // zeros in this case. See comment in Reset().
 221         return t.mat.Stride == 0
 222 }
 223
 224 // untranspose untransposes a matrix if applicable. If a is an Untransposer, then
 225 // untranspose returns the underlying matrix and true. If it is not, then it returns
 226 // the input matrix and false.
 227 func untransposeTri(a Triangular) (Triangular, bool) {
 228         if ut, ok := a.(UntransposeTrier); ok {
 229                 return ut.UntransposeTri(), true
 230         }
 231         return a, false
 232 }
 233
 234 // reuseAs resizes a zero receiver to an n×n triangular matrix with the given
 235 // orientation. If the receiver is non-zero, reuseAs checks that the receiver
 236 // is the correct size and orientation.
 237 func (t *TriDense) reuseAs(n int, kind TriKind) {
 238         ul := blas.Lower
 239         if kind == Upper {
 240                 ul = blas.Upper
 241         }
 242         if t.mat.N > t.cap {
 243                 panic(badTriCap)
 244         }
 245         if t.IsZero() {
 246                 t.mat = blas64.Triangular{
 247                         N:      n,
 248                         Stride: n,
 249                         Diag:   blas.NonUnit,
 250                         Data:   use(t.mat.Data, n*n),
 251                         Uplo:   ul,
 252                 }
 253                 t.cap = n
 254                 return
 255         }
 256         if t.mat.N != n {
 257                 panic(ErrShape)
 258         }
 259         if t.mat.Uplo != ul {
 260                 panic(ErrTriangle)
 261         }
 262 }
 263
 264 // isolatedWorkspace returns a new TriDense matrix w with the size of a and
 265 // returns a callback to defer which performs cleanup at the return of the call.
 266 // This should be used when a method receiver is the same pointer as an input argument.
 267 func (t *TriDense) isolatedWorkspace(a Triangular) (w *TriDense, restore func()) {
 268         n, kind := a.Triangle()
 269         if n == 0 {
 270                 panic("zero size")
 271         }
 272         w = getWorkspaceTri(n, kind, false)
 273         return w, func() {
 274                 t.Copy(w)
 275                 putWorkspaceTri(w)
 276         }
 277 }
 278
 279 // Copy makes a copy of elements of a into the receiver. It is similar to the
 280 // built-in copy; it copies as much as the overlap between the two matrices and
 281 // returns the number of rows and columns it copied. Only elements within the
 282 // receiver's non-zero triangle are set.
 283 //
 284 // See the Copier interface for more information.
 285 func (t *TriDense) Copy(a Matrix) (r, c int) {
 286         r, c = a.Dims()
 287         r = min(r, t.mat.N)
 288         c = min(c, t.mat.N)
 289         if r == 0 || c == 0 {
 290                 return 0, 0
 291         }
 292
 293         switch a := a.(type) {
 294         case RawMatrixer:
 295                 amat := a.RawMatrix()
 296                 if t.isUpper() {
 297                         for i := 0; i < r; i++ {
 298                                 copy(t.mat.Data[i*t.mat.Stride+i:i*t.mat.Stride+c], amat.Data[i*amat.Stride+i:i*amat.Stride+c])
 299                         }
 300                 } else {
 301                         for i := 0; i < r; i++ {
 302                                 copy(t.mat.Data[i*t.mat.Stride:i*t.mat.Stride+i+1], amat.Data[i*amat.Stride:i*amat.Stride+i+1])
 303                         }
 304                 }
 305         case RawTriangular:
 306                 amat := a.RawTriangular()
 307                 aIsUpper := isUpperUplo(amat.Uplo)
 308                 tIsUpper := t.isUpper()
 309                 switch {
 310                 case tIsUpper && aIsUpper:
 311                         for i := 0; i < r; i++ {
 312                                 copy(t.mat.Data[i*t.mat.Stride+i:i*t.mat.Stride+c], amat.Data[i*amat.Stride+i:i*amat.Stride+c])
 313                         }
 314                 case !tIsUpper && !aIsUpper:
 315                         for i := 0; i < r; i++ {
 316                                 copy(t.mat.Data[i*t.mat.Stride:i*t.mat.Stride+i+1], amat.Data[i*amat.Stride:i*amat.Stride+i+1])
 317                         }
 318                 default:
 319                         for i := 0; i < r; i++ {
 320                                 t.set(i, i, amat.Data[i*amat.Stride+i])
 321                         }
 322                 }
 323         default:
 324                 isUpper := t.isUpper()
 325                 for i := 0; i < r; i++ {
 326                         if isUpper {
 327                                 for j := i; j < c; j++ {
 328                                         t.set(i, j, a.At(i, j))
 329                                 }
 330                         } else {
 331                                 for j := 0; j <= i; j++ {
 332                                         t.set(i, j, a.At(i, j))
 333                                 }
 334                         }
 335                 }
 336         }
 337
 338         return r, c
 339 }
 340
 341 // InverseTri computes the inverse of the triangular matrix a, storing the result
 342 // into the receiver. If a is ill-conditioned, a Condition error will be returned.
 343 // Note that matrix inversion is numerically unstable, and should generally be
 344 // avoided where possible, for example by using the Solve routines.
 345 func (t *TriDense) InverseTri(a Triangular) error {
 346         if rt, ok := a.(RawTriangular); ok {
 347                 t.checkOverlap(generalFromTriangular(rt.RawTriangular()))
 348         }
 349         n, _ := a.Triangle()
 350         t.reuseAs(a.Triangle())
 351         t.Copy(a)
 352         work := getFloats(3*n, false)
 353         iwork := getInts(n, false)
 354         cond := lapack64.Trcon(CondNorm, t.mat, work, iwork)
 355         putFloats(work)
 356         putInts(iwork)
 357         if math.IsInf(cond, 1) {
 358                 return Condition(cond)
 359         }
 360         ok := lapack64.Trtri(t.mat)
 361         if !ok {
 362                 return Condition(math.Inf(1))
 363         }
 364         if cond > ConditionTolerance {
 365                 return Condition(cond)
 366         }
 367         return nil
 368 }
 369
 370 // MulTri takes the product of triangular matrices a and b and places the result
 371 // in the receiver. The size of a and b must match, and they both must have the
 372 // same TriKind, or Mul will panic.
 373 func (t *TriDense) MulTri(a, b Triangular) {
 374         n, kind := a.Triangle()
 375         nb, kindb := b.Triangle()
 376         if n != nb {
 377                 panic(ErrShape)
 378         }
 379         if kind != kindb {
 380                 panic(ErrTriangle)
 381         }
 382
 383         aU, _ := untransposeTri(a)
 384         bU, _ := untransposeTri(b)
 385         t.reuseAs(n, kind)
 386         var restore func()
 387         if t == aU {
 388                 t, restore = t.isolatedWorkspace(aU)
 389                 defer restore()
 390         } else if t == bU {
 391                 t, restore = t.isolatedWorkspace(bU)
 392                 defer restore()
 393         }
 394
 395         // TODO(btracey): Improve the set of fast-paths.
 396         if kind == Upper {
 397                 for i := 0; i < n; i++ {
 398                         for j := i; j < n; j++ {
 399                                 var v float64
 400                                 for k := i; k <= j; k++ {
 401                                         v += a.At(i, k) * b.At(k, j)
 402                                 }
 403                                 t.SetTri(i, j, v)
 404                         }
 405                 }
 406                 return
 407         }
 408         for i := 0; i < n; i++ {
 409                 for j := 0; j <= i; j++ {
 410                         var v float64
 411                         for k := j; k <= i; k++ {
 412                                 v += a.At(i, k) * b.At(k, j)
 413                         }
 414                         t.SetTri(i, j, v)
 415                 }
 416         }
 417 }
 418
 419 // ScaleTri multiplies the elements of a by f, placing the result in the receiver.
 420 // If the receiver is non-zero, the size and kind of the receiver must match
 421 // the input, or ScaleTri will panic.
 422 func (t *TriDense) ScaleTri(f float64, a Triangular) {
 423         n, kind := a.Triangle()
 424         t.reuseAs(n, kind)
 425
 426         // TODO(btracey): Improve the set of fast-paths.
 427         switch a := a.(type) {
 428         case RawTriangular:
 429                 amat := a.RawTriangular()
 430                 if kind == Upper {
 431                         for i := 0; i < n; i++ {
 432                                 ts := t.mat.Data[i*t.mat.Stride+i : i*t.mat.Stride+n]
 433                                 as := amat.Data[i*amat.Stride+i : i*amat.Stride+n]
 434                                 for i, v := range as {
 435                                         ts[i] = v * f
 436                                 }
 437                         }
 438                         return
 439                 }
 440                 for i := 0; i < n; i++ {
 441                         ts := t.mat.Data[i*t.mat.Stride : i*t.mat.Stride+i+1]
 442                         as := amat.Data[i*amat.Stride : i*amat.Stride+i+1]
 443                         for i, v := range as {
 444                                 ts[i] = v * f
 445                         }
 446                 }
 447                 return
 448         default:
 449                 isUpper := kind == Upper
 450                 for i := 0; i < n; i++ {
 451                         if isUpper {
 452                                 for j := i; j < n; j++ {
 453                                         t.set(i, j, f*a.At(i, j))
 454                                 }
 455                         } else {
 456                                 for j := 0; j <= i; j++ {
 457                                         t.set(i, j, f*a.At(i, j))
 458                                 }
 459                         }
 460                 }
 461         }
 462 }
 463
 464 // copySymIntoTriangle copies a symmetric matrix into a TriDense
 465 func copySymIntoTriangle(t *TriDense, s Symmetric) {
 466         n, upper := t.Triangle()
 467         ns := s.Symmetric()
 468         if n != ns {
 469                 panic("mat: triangle size mismatch")
 470         }
 471         ts := t.mat.Stride
 472         if rs, ok := s.(RawSymmetricer); ok {
 473                 sd := rs.RawSymmetric()
 474                 ss := sd.Stride
 475                 if upper {
 476                         if sd.Uplo == blas.Upper {
 477                                 for i := 0; i < n; i++ {
 478                                         copy(t.mat.Data[i*ts+i:i*ts+n], sd.Data[i*ss+i:i*ss+n])
 479                                 }
 480                                 return
 481                         }
 482                         for i := 0; i < n; i++ {
 483                                 for j := i; j < n; j++ {
 484                                         t.mat.Data[i*ts+j] = sd.Data[j*ss+i]
 485                                 }
 486                         }
 487                         return
 488                 }
 489                 if sd.Uplo == blas.Upper {
 490                         for i := 0; i < n; i++ {
 491                                 for j := 0; j <= i; j++ {
 492                                         t.mat.Data[i*ts+j] = sd.Data[j*ss+i]
 493                                 }
 494                         }
 495                         return
 496                 }
 497                 for i := 0; i < n; i++ {
 498                         copy(t.mat.Data[i*ts:i*ts+i+1], sd.Data[i*ss:i*ss+i+1])
 499                 }
 500                 return
 501         }
 502         if upper {
 503                 for i := 0; i < n; i++ {
 504                         for j := i; j < n; j++ {
 505                                 t.mat.Data[i*ts+j] = s.At(i, j)
 506                         }
 507                 }
 508                 return
 509         }
 510         for i := 0; i < n; i++ {
 511                 for j := 0; j <= i; j++ {
 512                         t.mat.Data[i*ts+j] = s.At(i, j)
 513                 }
 514         }
 515 }
 516
 517 // DoNonZero calls the function fn for each of the non-zero elements of t. The function fn
 518 // takes a row/column index and the element value of t at (i, j).
 519 func (t *TriDense) DoNonZero(fn func(i, j int, v float64)) {
 520         if t.isUpper() {
 521                 for i := 0; i < t.mat.N; i++ {
 522                         for j := i; j < t.mat.N; j++ {
 523                                 v := t.at(i, j)
 524                                 if v != 0 {
 525                                         fn(i, j, v)
 526                                 }
 527                         }
 528                 }
 529                 return
 530         }
 531         for i := 0; i < t.mat.N; i++ {
 532                 for j := 0; j <= i; j++ {
 533                         v := t.at(i, j)
 534                         if v != 0 {
 535                                 fn(i, j, v)
 536                         }
 537                 }
 538         }
 539 }
 540
 541 // DoRowNonZero calls the function fn for each of the non-zero elements of row i of t. The function fn
 542 // takes a row/column index and the element value of t at (i, j).
 543 func (t *TriDense) DoRowNonZero(i int, fn func(i, j int, v float64)) {
 544         if i < 0 || t.mat.N <= i {
 545                 panic(ErrRowAccess)
 546         }
 547         if t.isUpper() {
 548                 for j := i; j < t.mat.N; j++ {
 549                         v := t.at(i, j)
 550                         if v != 0 {
 551                                 fn(i, j, v)
 552                         }
 553                 }
 554                 return
 555         }
 556         for j := 0; j <= i; j++ {
 557                 v := t.at(i, j)
 558                 if v != 0 {
 559                         fn(i, j, v)
 560                 }
 561         }
 562 }
 563
 564 // DoColNonZero calls the function fn for each of the non-zero elements of column j of t. The function fn
 565 // takes a row/column index and the element value of t at (i, j).
 566 func (t *TriDense) DoColNonZero(j int, fn func(i, j int, v float64)) {
 567         if j < 0 || t.mat.N <= j {
 568                 panic(ErrColAccess)
 569         }
 570         if t.isUpper() {
 571                 for i := 0; i <= j; i++ {
 572                         v := t.at(i, j)
 573                         if v != 0 {
 574                                 fn(i, j, v)
 575                         }
 576                 }
 577                 return
 578         }
 579         for i := j; i < t.mat.N; i++ {
 580                 v := t.at(i, j)
 581                 if v != 0 {
 582                         fn(i, j, v)
 583                 }
 584         }
 585 }