vendor/gonum.org/v1/gonum/mat/symmetric_test.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         "fmt"
   9         "os"
  10         "reflect"
  11         "testing"
  12
  13         "golang.org/x/exp/rand"
  14
  15         "gonum.org/v1/gonum/blas"
  16         "gonum.org/v1/gonum/blas/blas64"
  17         "gonum.org/v1/gonum/floats"
  18 )
  19
  20 func TestNewSymmetric(t *testing.T) {
  21         for i, test := range []struct {
  22                 data []float64
  23                 n    int
  24                 mat  *SymDense
  25         }{
  26                 {
  27                         data: []float64{
  28                                 1, 2, 3,
  29                                 4, 5, 6,
  30                                 7, 8, 9,
  31                         },
  32                         n: 3,
  33                         mat: &SymDense{
  34                                 mat: blas64.Symmetric{
  35                                         N:      3,
  36                                         Stride: 3,
  37                                         Uplo:   blas.Upper,
  38                                         Data:   []float64{1, 2, 3, 4, 5, 6, 7, 8, 9},
  39                                 },
  40                                 cap: 3,
  41                         },
  42                 },
  43         } {
  44                 sym := NewSymDense(test.n, test.data)
  45                 rows, cols := sym.Dims()
  46
  47                 if rows != test.n {
  48                         t.Errorf("unexpected number of rows for test %d: got: %d want: %d", i, rows, test.n)
  49                 }
  50                 if cols != test.n {
  51                         t.Errorf("unexpected number of cols for test %d: got: %d want: %d", i, cols, test.n)
  52                 }
  53                 if !reflect.DeepEqual(sym, test.mat) {
  54                         t.Errorf("unexpected data slice for test %d: got: %v want: %v", i, sym, test.mat)
  55                 }
  56
  57                 m := NewDense(test.n, test.n, test.data)
  58                 if !reflect.DeepEqual(sym.mat.Data, m.mat.Data) {
  59                         t.Errorf("unexpected data slice mismatch for test %d: got: %v want: %v", i, sym.mat.Data, m.mat.Data)
  60                 }
  61         }
  62
  63         panicked, message := panics(func() { NewSymDense(3, []float64{1, 2}) })
  64         if !panicked || message != ErrShape.Error() {
  65                 t.Error("expected panic for invalid data slice length")
  66         }
  67 }
  68
  69 func TestSymAtSet(t *testing.T) {
  70         sym := &SymDense{
  71                 mat: blas64.Symmetric{
  72                         N:      3,
  73                         Stride: 3,
  74                         Uplo:   blas.Upper,
  75                         Data:   []float64{1, 2, 3, 4, 5, 6, 7, 8, 9},
  76                 },
  77                 cap: 3,
  78         }
  79         rows, cols := sym.Dims()
  80
  81         // Check At out of bounds
  82         for _, row := range []int{-1, rows, rows + 1} {
  83                 panicked, message := panics(func() { sym.At(row, 0) })
  84                 if !panicked || message != ErrRowAccess.Error() {
  85                         t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
  86                 }
  87         }
  88         for _, col := range []int{-1, cols, cols + 1} {
  89                 panicked, message := panics(func() { sym.At(0, col) })
  90                 if !panicked || message != ErrColAccess.Error() {
  91                         t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
  92                 }
  93         }
  94
  95         // Check Set out of bounds
  96         for _, row := range []int{-1, rows, rows + 1} {
  97                 panicked, message := panics(func() { sym.SetSym(row, 0, 1.2) })
  98                 if !panicked || message != ErrRowAccess.Error() {
  99                         t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
 100                 }
 101         }
 102         for _, col := range []int{-1, cols, cols + 1} {
 103                 panicked, message := panics(func() { sym.SetSym(0, col, 1.2) })
 104                 if !panicked || message != ErrColAccess.Error() {
 105                         t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
 106                 }
 107         }
 108
 109         for _, st := range []struct {
 110                 row, col  int
 111                 orig, new float64
 112         }{
 113                 {row: 1, col: 2, orig: 6, new: 15},
 114                 {row: 2, col: 1, orig: 15, new: 12},
 115         } {
 116                 if e := sym.At(st.row, st.col); e != st.orig {
 117                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.row, st.col, e, st.orig)
 118                 }
 119                 if e := sym.At(st.col, st.row); e != st.orig {
 120                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.col, st.row, e, st.orig)
 121                 }
 122                 sym.SetSym(st.row, st.col, st.new)
 123                 if e := sym.At(st.row, st.col); e != st.new {
 124                         t.Errorf("unexpected value for At(%d, %d) after SetSym(%[1]d, %[2]d, %[4]v): got: %[3]v want: %v", st.row, st.col, e, st.new)
 125                 }
 126                 if e := sym.At(st.col, st.row); e != st.new {
 127                         t.Errorf("unexpected value for At(%d, %d) after SetSym(%[2]d, %[1]d, %[4]v): got: %[3]v want: %v", st.col, st.row, e, st.new)
 128                 }
 129         }
 130 }
 131
 132 func TestSymAdd(t *testing.T) {
 133         for _, test := range []struct {
 134                 n int
 135         }{
 136                 {n: 1},
 137                 {n: 2},
 138                 {n: 3},
 139                 {n: 4},
 140                 {n: 5},
 141                 {n: 10},
 142         } {
 143                 n := test.n
 144                 a := NewSymDense(n, nil)
 145                 for i := range a.mat.Data {
 146                         a.mat.Data[i] = rand.Float64()
 147                 }
 148                 b := NewSymDense(n, nil)
 149                 for i := range a.mat.Data {
 150                         b.mat.Data[i] = rand.Float64()
 151                 }
 152                 var m Dense
 153                 m.Add(a, b)
 154
 155                 // Check with new receiver
 156                 var s SymDense
 157                 s.AddSym(a, b)
 158                 for i := 0; i < n; i++ {
 159                         for j := i; j < n; j++ {
 160                                 want := m.At(i, j)
 161                                 if got := s.At(i, j); got != want {
 162                                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
 163                                 }
 164                         }
 165                 }
 166
 167                 // Check with equal receiver
 168                 s.CopySym(a)
 169                 s.AddSym(&s, b)
 170                 for i := 0; i < n; i++ {
 171                         for j := i; j < n; j++ {
 172                                 want := m.At(i, j)
 173                                 if got := s.At(i, j); got != want {
 174                                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
 175                                 }
 176                         }
 177                 }
 178         }
 179
 180         method := func(receiver, a, b Matrix) {
 181                 type addSymer interface {
 182                         AddSym(a, b Symmetric)
 183                 }
 184                 rd := receiver.(addSymer)
 185                 rd.AddSym(a.(Symmetric), b.(Symmetric))
 186         }
 187         denseComparison := func(receiver, a, b *Dense) {
 188                 receiver.Add(a, b)
 189         }
 190         testTwoInput(t, "AddSym", &SymDense{}, method, denseComparison, legalTypesSym, legalSizeSameSquare, 1e-14)
 191 }
 192
 193 func TestCopy(t *testing.T) {
 194         for _, test := range []struct {
 195                 n int
 196         }{
 197                 {n: 1},
 198                 {n: 2},
 199                 {n: 3},
 200                 {n: 4},
 201                 {n: 5},
 202                 {n: 10},
 203         } {
 204                 n := test.n
 205                 a := NewSymDense(n, nil)
 206                 for i := range a.mat.Data {
 207                         a.mat.Data[i] = rand.Float64()
 208                 }
 209                 s := NewSymDense(n, nil)
 210                 s.CopySym(a)
 211                 for i := 0; i < n; i++ {
 212                         for j := i; j < n; j++ {
 213                                 want := a.At(i, j)
 214                                 if got := s.At(i, j); got != want {
 215                                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
 216                                 }
 217                         }
 218                 }
 219         }
 220 }
 221
 222 // TODO(kortschak) Roll this into testOneInput when it exists.
 223 // https://github.com/gonum/matrix/issues/171
 224 func TestSymCopyPanic(t *testing.T) {
 225         var (
 226                 a SymDense
 227                 n int
 228         )
 229         m := NewSymDense(1, nil)
 230         panicked, message := panics(func() { n = m.CopySym(&a) })
 231         if panicked {
 232                 t.Errorf("unexpected panic: %v", message)
 233         }
 234         if n != 0 {
 235                 t.Errorf("unexpected n: got: %d want: 0", n)
 236         }
 237 }
 238
 239 func TestSymRankOne(t *testing.T) {
 240         for _, test := range []struct {
 241                 n int
 242         }{
 243                 {n: 1},
 244                 {n: 2},
 245                 {n: 3},
 246                 {n: 4},
 247                 {n: 5},
 248                 {n: 10},
 249         } {
 250                 n := test.n
 251                 alpha := 2.0
 252                 a := NewSymDense(n, nil)
 253                 for i := range a.mat.Data {
 254                         a.mat.Data[i] = rand.Float64()
 255                 }
 256                 x := make([]float64, n)
 257                 for i := range x {
 258                         x[i] = rand.Float64()
 259                 }
 260
 261                 xMat := NewDense(n, 1, x)
 262                 var m Dense
 263                 m.Mul(xMat, xMat.T())
 264                 m.Scale(alpha, &m)
 265                 m.Add(&m, a)
 266
 267                 // Check with new receiver
 268                 s := NewSymDense(n, nil)
 269                 s.SymRankOne(a, alpha, NewVecDense(len(x), x))
 270                 for i := 0; i < n; i++ {
 271                         for j := i; j < n; j++ {
 272                                 want := m.At(i, j)
 273                                 if got := s.At(i, j); got != want {
 274                                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
 275                                 }
 276                         }
 277                 }
 278
 279                 // Check with reused receiver
 280                 copy(s.mat.Data, a.mat.Data)
 281                 s.SymRankOne(s, alpha, NewVecDense(len(x), x))
 282                 for i := 0; i < n; i++ {
 283                         for j := i; j < n; j++ {
 284                                 want := m.At(i, j)
 285                                 if got := s.At(i, j); got != want {
 286                                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
 287                                 }
 288                         }
 289                 }
 290         }
 291
 292         alpha := 3.0
 293         method := func(receiver, a, b Matrix) {
 294                 type SymRankOner interface {
 295                         SymRankOne(a Symmetric, alpha float64, x Vector)
 296                 }
 297                 rd := receiver.(SymRankOner)
 298                 rd.SymRankOne(a.(Symmetric), alpha, b.(Vector))
 299         }
 300         denseComparison := func(receiver, a, b *Dense) {
 301                 var tmp Dense
 302                 tmp.Mul(b, b.T())
 303                 tmp.Scale(alpha, &tmp)
 304                 receiver.Add(a, &tmp)
 305         }
 306         legalTypes := func(a, b Matrix) bool {
 307                 _, ok := a.(Symmetric)
 308                 if !ok {
 309                         return false
 310                 }
 311                 _, ok = b.(Vector)
 312                 return ok
 313         }
 314         legalSize := func(ar, ac, br, bc int) bool {
 315                 if ar != ac {
 316                         return false
 317                 }
 318                 return br == ar
 319         }
 320         testTwoInput(t, "SymRankOne", &SymDense{}, method, denseComparison, legalTypes, legalSize, 1e-14)
 321 }
 322
 323 func TestIssue250SymRankOne(t *testing.T) {
 324         x := NewVecDense(5, []float64{1, 2, 3, 4, 5})
 325         var s1, s2 SymDense
 326         s1.SymRankOne(NewSymDense(5, nil), 1, x)
 327         s2.SymRankOne(NewSymDense(5, nil), 1, x)
 328         s2.SymRankOne(NewSymDense(5, nil), 1, x)
 329         if !Equal(&s1, &s2) {
 330                 t.Error("unexpected result from repeat")
 331         }
 332 }
 333
 334 func TestRankTwo(t *testing.T) {
 335         for _, test := range []struct {
 336                 n int
 337         }{
 338                 {n: 1},
 339                 {n: 2},
 340                 {n: 3},
 341                 {n: 4},
 342                 {n: 5},
 343                 {n: 10},
 344         } {
 345                 n := test.n
 346                 alpha := 2.0
 347                 a := NewSymDense(n, nil)
 348                 for i := range a.mat.Data {
 349                         a.mat.Data[i] = rand.Float64()
 350                 }
 351                 x := make([]float64, n)
 352                 y := make([]float64, n)
 353                 for i := range x {
 354                         x[i] = rand.Float64()
 355                         y[i] = rand.Float64()
 356                 }
 357
 358                 xMat := NewDense(n, 1, x)
 359                 yMat := NewDense(n, 1, y)
 360                 var m Dense
 361                 m.Mul(xMat, yMat.T())
 362                 var tmp Dense
 363                 tmp.Mul(yMat, xMat.T())
 364                 m.Add(&m, &tmp)
 365                 m.Scale(alpha, &m)
 366                 m.Add(&m, a)
 367
 368                 // Check with new receiver
 369                 s := NewSymDense(n, nil)
 370                 s.RankTwo(a, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y))
 371                 for i := 0; i < n; i++ {
 372                         for j := i; j < n; j++ {
 373                                 if !floats.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
 374                                         t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j))
 375                                 }
 376                         }
 377                 }
 378
 379                 // Check with reused receiver
 380                 copy(s.mat.Data, a.mat.Data)
 381                 s.RankTwo(s, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y))
 382                 for i := 0; i < n; i++ {
 383                         for j := i; j < n; j++ {
 384                                 if !floats.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
 385                                         t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j))
 386                                 }
 387                         }
 388                 }
 389         }
 390 }
 391
 392 func TestSymRankK(t *testing.T) {
 393         alpha := 3.0
 394         method := func(receiver, a, b Matrix) {
 395                 type SymRankKer interface {
 396                         SymRankK(a Symmetric, alpha float64, x Matrix)
 397                 }
 398                 rd := receiver.(SymRankKer)
 399                 rd.SymRankK(a.(Symmetric), alpha, b)
 400         }
 401         denseComparison := func(receiver, a, b *Dense) {
 402                 var tmp Dense
 403                 tmp.Mul(b, b.T())
 404                 tmp.Scale(alpha, &tmp)
 405                 receiver.Add(a, &tmp)
 406         }
 407         legalTypes := func(a, b Matrix) bool {
 408                 _, ok := a.(Symmetric)
 409                 return ok
 410         }
 411         legalSize := func(ar, ac, br, bc int) bool {
 412                 if ar != ac {
 413                         return false
 414                 }
 415                 return br == ar
 416         }
 417         testTwoInput(t, "SymRankK", &SymDense{}, method, denseComparison, legalTypes, legalSize, 1e-14)
 418 }
 419
 420 func TestSymOuterK(t *testing.T) {
 421         for _, f := range []float64{0.5, 1, 3} {
 422                 method := func(receiver, x Matrix) {
 423                         type SymOuterKer interface {
 424                                 SymOuterK(alpha float64, x Matrix)
 425                         }
 426                         rd := receiver.(SymOuterKer)
 427                         rd.SymOuterK(f, x)
 428                 }
 429                 denseComparison := func(receiver, x *Dense) {
 430                         receiver.Mul(x, x.T())
 431                         receiver.Scale(f, receiver)
 432                 }
 433                 testOneInput(t, "SymOuterK", &SymDense{}, method, denseComparison, isAnyType, isAnySize, 1e-14)
 434         }
 435 }
 436
 437 func TestIssue250SymOuterK(t *testing.T) {
 438         x := NewVecDense(5, []float64{1, 2, 3, 4, 5})
 439         var s1, s2 SymDense
 440         s1.SymOuterK(1, x)
 441         s2.SymOuterK(1, x)
 442         s2.SymOuterK(1, x)
 443         if !Equal(&s1, &s2) {
 444                 t.Error("unexpected result from repeat")
 445         }
 446 }
 447
 448 func TestScaleSym(t *testing.T) {
 449         for _, f := range []float64{0.5, 1, 3} {
 450                 method := func(receiver, a Matrix) {
 451                         type ScaleSymer interface {
 452                                 ScaleSym(f float64, a Symmetric)
 453                         }
 454                         rd := receiver.(ScaleSymer)
 455                         rd.ScaleSym(f, a.(Symmetric))
 456                 }
 457                 denseComparison := func(receiver, a *Dense) {
 458                         receiver.Scale(f, a)
 459                 }
 460                 testOneInput(t, "ScaleSym", &SymDense{}, method, denseComparison, legalTypeSym, isSquare, 1e-14)
 461         }
 462 }
 463
 464 func TestSubsetSym(t *testing.T) {
 465         for _, test := range []struct {
 466                 a    *SymDense
 467                 dims []int
 468                 ans  *SymDense
 469         }{
 470                 {
 471                         a: NewSymDense(3, []float64{
 472                                 1, 2, 3,
 473                                 0, 4, 5,
 474                                 0, 0, 6,
 475                         }),
 476                         dims: []int{0, 2},
 477                         ans: NewSymDense(2, []float64{
 478                                 1, 3,
 479                                 0, 6,
 480                         }),
 481                 },
 482                 {
 483                         a: NewSymDense(3, []float64{
 484                                 1, 2, 3,
 485                                 0, 4, 5,
 486                                 0, 0, 6,
 487                         }),
 488                         dims: []int{2, 0},
 489                         ans: NewSymDense(2, []float64{
 490                                 6, 3,
 491                                 0, 1,
 492                         }),
 493                 },
 494                 {
 495                         a: NewSymDense(3, []float64{
 496                                 1, 2, 3,
 497                                 0, 4, 5,
 498                                 0, 0, 6,
 499                         }),
 500                         dims: []int{1, 1, 1},
 501                         ans: NewSymDense(3, []float64{
 502                                 4, 4, 4,
 503                                 0, 4, 4,
 504                                 0, 0, 4,
 505                         }),
 506                 },
 507         } {
 508                 var s SymDense
 509                 s.SubsetSym(test.a, test.dims)
 510                 if !Equal(&s, test.ans) {
 511                         t.Errorf("SubsetSym mismatch dims %v\nGot:\n% v\nWant:\n% v\n", test.dims, s, test.ans)
 512                 }
 513         }
 514
 515         dims := []int{0, 2}
 516         maxDim := dims[0]
 517         for _, v := range dims {
 518                 if maxDim < v {
 519                         maxDim = v
 520                 }
 521         }
 522         method := func(receiver, a Matrix) {
 523                 type SubsetSymer interface {
 524                         SubsetSym(a Symmetric, set []int)
 525                 }
 526                 rd := receiver.(SubsetSymer)
 527                 rd.SubsetSym(a.(Symmetric), dims)
 528         }
 529         denseComparison := func(receiver, a *Dense) {
 530                 *receiver = *NewDense(len(dims), len(dims), nil)
 531                 sz := len(dims)
 532                 for i := 0; i < sz; i++ {
 533                         for j := 0; j < sz; j++ {
 534                                 receiver.Set(i, j, a.At(dims[i], dims[j]))
 535                         }
 536                 }
 537         }
 538         legalSize := func(ar, ac int) bool {
 539                 return ar == ac && ar > maxDim
 540         }
 541
 542         testOneInput(t, "SubsetSym", &SymDense{}, method, denseComparison, legalTypeSym, legalSize, 0)
 543 }
 544
 545 func TestViewGrowSquare(t *testing.T) {
 546         // n is the size of the original SymDense.
 547         // The first view uses start1, span1. The second view uses start2, span2 on
 548         // the first view.
 549         for _, test := range []struct {
 550                 n, start1, span1, start2, span2 int
 551         }{
 552                 {10, 0, 10, 0, 10},
 553                 {10, 0, 8, 0, 8},
 554                 {10, 2, 8, 0, 6},
 555                 {10, 2, 7, 4, 2},
 556                 {10, 2, 6, 0, 5},
 557         } {
 558                 n := test.n
 559                 s := NewSymDense(n, nil)
 560                 for i := 0; i < n; i++ {
 561                         for j := i; j < n; j++ {
 562                                 s.SetSym(i, j, float64((i+1)*n+j+1))
 563                         }
 564                 }
 565
 566                 // Take a subset and check the view matches.
 567                 start1 := test.start1
 568                 span1 := test.span1
 569                 v := s.SliceSquare(start1, start1+span1).(*SymDense)
 570                 for i := 0; i < span1; i++ {
 571                         for j := i; j < span1; j++ {
 572                                 if v.At(i, j) != s.At(start1+i, start1+j) {
 573                                         t.Errorf("View mismatch")
 574                                 }
 575                         }
 576                 }
 577
 578                 start2 := test.start2
 579                 span2 := test.span2
 580                 v2 := v.SliceSquare(start2, start2+span2).(*SymDense)
 581
 582                 for i := 0; i < span2; i++ {
 583                         for j := i; j < span2; j++ {
 584                                 if v2.At(i, j) != s.At(start1+start2+i, start1+start2+j) {
 585                                         t.Errorf("Second view mismatch")
 586                                 }
 587                         }
 588                 }
 589
 590                 // Check that a write to the view is reflected in the original.
 591                 v2.SetSym(0, 0, 1.2)
 592                 if s.At(start1+start2, start1+start2) != 1.2 {
 593                         t.Errorf("Write to view not reflected in original")
 594                 }
 595
 596                 // Grow the matrix back to the original view
 597                 gn := n - start1 - start2
 598                 g := v2.GrowSquare(gn - v2.Symmetric()).(*SymDense)
 599                 g.SetSym(1, 1, 2.2)
 600
 601                 for i := 0; i < gn; i++ {
 602                         for j := 0; j < gn; j++ {
 603                                 if g.At(i, j) != s.At(start1+start2+i, start1+start2+j) {
 604                                         t.Errorf("Grow mismatch")
 605
 606                                         fmt.Printf("g=\n% v\n", Formatted(g))
 607                                         fmt.Printf("s=\n% v\n", Formatted(s))
 608                                         os.Exit(1)
 609                                 }
 610                         }
 611                 }
 612
 613                 // View g, then grow it and make sure all the elements were copied.
 614                 gv := g.SliceSquare(0, gn-1).(*SymDense)
 615
 616                 gg := gv.GrowSquare(2)
 617                 for i := 0; i < gn; i++ {
 618                         for j := 0; j < gn; j++ {
 619                                 if g.At(i, j) != gg.At(i, j) {
 620                                         t.Errorf("Expand mismatch")
 621                                 }
 622                         }
 623                 }
 624         }
 625 }
 626
 627 func TestPowPSD(t *testing.T) {
 628         for cas, test := range []struct {
 629                 a   *SymDense
 630                 pow float64
 631                 ans *SymDense
 632         }{
 633                 // Comparison with Matlab.
 634                 {
 635                         a:   NewSymDense(2, []float64{10, 5, 5, 12}),
 636                         pow: 0.5,
 637                         ans: NewSymDense(2, []float64{3.065533767740645, 0.776210486171016, 0.776210486171016, 3.376017962209052}),
 638                 },
 639                 {
 640                         a:   NewSymDense(2, []float64{11, -1, -1, 8}),
 641                         pow: 0.5,
 642                         ans: NewSymDense(2, []float64{3.312618742210524, -0.162963396980939, -0.162963396980939, 2.823728551267709}),
 643                 },
 644                 {
 645                         a:   NewSymDense(2, []float64{10, 5, 5, 12}),
 646                         pow: -0.5,
 647                         ans: NewSymDense(2, []float64{0.346372134547712, -0.079637515547296, -0.079637515547296, 0.314517128328794}),
 648                 },
 649                 {
 650                         a:   NewSymDense(3, []float64{15, -1, -3, -1, 8, 6, -3, 6, 14}),
 651                         pow: 0.6,
 652                         ans: NewSymDense(3, []float64{
 653                                 5.051214323034288, -0.163162161893975, -0.612153996497505,
 654                                 -0.163162161893976, 3.283474884617009, 1.432842761381493,
 655                                 -0.612153996497505, 1.432842761381494, 4.695873060862573,
 656                         }),
 657                 },
 658         } {
 659                 var s SymDense
 660                 err := s.PowPSD(test.a, test.pow)
 661                 if err != nil {
 662                         panic("bad test")
 663                 }
 664                 if !EqualApprox(&s, test.ans, 1e-10) {
 665                         t.Errorf("Case %d, pow mismatch", cas)
 666                         fmt.Println(Formatted(&s))
 667                         fmt.Println(Formatted(test.ans))
 668                 }
 669         }
 670
 671         // Compare with Dense.Pow
 672         rnd := rand.New(rand.NewSource(1))
 673         for dim := 2; dim < 10; dim++ {
 674                 for pow := 2; pow < 6; pow++ {
 675                         a := NewDense(dim, dim, nil)
 676                         for i := 0; i < dim; i++ {
 677                                 for j := 0; j < dim; j++ {
 678                                         a.Set(i, j, rnd.Float64())
 679                                 }
 680                         }
 681                         var mat SymDense
 682                         mat.SymOuterK(1, a)
 683
 684                         var sym SymDense
 685                         sym.PowPSD(&mat, float64(pow))
 686
 687                         var dense Dense
 688                         dense.Pow(&mat, pow)
 689
 690                         if !EqualApprox(&sym, &dense, 1e-10) {
 691                                 t.Errorf("Dim %d: pow mismatch", dim)
 692                         }
 693                 }
 694         }
 695 }