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

   1 // Copyright ©2013 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         "math"
  10         "reflect"
  11         "testing"
  12
  13         "golang.org/x/exp/rand"
  14
  15         "gonum.org/v1/gonum/blas/blas64"
  16         "gonum.org/v1/gonum/floats"
  17 )
  18
  19 func asBasicMatrix(d *Dense) Matrix            { return (*basicMatrix)(d) }
  20 func asBasicSymmetric(s *SymDense) Matrix      { return (*basicSymmetric)(s) }
  21 func asBasicTriangular(t *TriDense) Triangular { return (*basicTriangular)(t) }
  22
  23 func TestNewDense(t *testing.T) {
  24         for i, test := range []struct {
  25                 a          []float64
  26                 rows, cols int
  27                 min, max   float64
  28                 fro        float64
  29                 mat        *Dense
  30         }{
  31                 {
  32                         []float64{
  33                                 0, 0, 0,
  34                                 0, 0, 0,
  35                                 0, 0, 0,
  36                         },
  37                         3, 3,
  38                         0, 0,
  39                         0,
  40                         &Dense{
  41                                 mat: blas64.General{
  42                                         Rows: 3, Cols: 3,
  43                                         Stride: 3,
  44                                         Data:   []float64{0, 0, 0, 0, 0, 0, 0, 0, 0},
  45                                 },
  46                                 capRows: 3, capCols: 3,
  47                         },
  48                 },
  49                 {
  50                         []float64{
  51                                 1, 1, 1,
  52                                 1, 1, 1,
  53                                 1, 1, 1,
  54                         },
  55                         3, 3,
  56                         1, 1,
  57                         3,
  58                         &Dense{
  59                                 mat: blas64.General{
  60                                         Rows: 3, Cols: 3,
  61                                         Stride: 3,
  62                                         Data:   []float64{1, 1, 1, 1, 1, 1, 1, 1, 1},
  63                                 },
  64                                 capRows: 3, capCols: 3,
  65                         },
  66                 },
  67                 {
  68                         []float64{
  69                                 1, 0, 0,
  70                                 0, 1, 0,
  71                                 0, 0, 1,
  72                         },
  73                         3, 3,
  74                         0, 1,
  75                         1.7320508075688772,
  76                         &Dense{
  77                                 mat: blas64.General{
  78                                         Rows: 3, Cols: 3,
  79                                         Stride: 3,
  80                                         Data:   []float64{1, 0, 0, 0, 1, 0, 0, 0, 1},
  81                                 },
  82                                 capRows: 3, capCols: 3,
  83                         },
  84                 },
  85                 {
  86                         []float64{
  87                                 -1, 0, 0,
  88                                 0, -1, 0,
  89                                 0, 0, -1,
  90                         },
  91                         3, 3,
  92                         -1, 0,
  93                         1.7320508075688772,
  94                         &Dense{
  95                                 mat: blas64.General{
  96                                         Rows: 3, Cols: 3,
  97                                         Stride: 3,
  98                                         Data:   []float64{-1, 0, 0, 0, -1, 0, 0, 0, -1},
  99                                 },
 100                                 capRows: 3, capCols: 3,
 101                         },
 102                 },
 103                 {
 104                         []float64{
 105                                 1, 2, 3,
 106                                 4, 5, 6,
 107                         },
 108                         2, 3,
 109                         1, 6,
 110                         9.539392014169458,
 111                         &Dense{
 112                                 mat: blas64.General{
 113                                         Rows: 2, Cols: 3,
 114                                         Stride: 3,
 115                                         Data:   []float64{1, 2, 3, 4, 5, 6},
 116                                 },
 117                                 capRows: 2, capCols: 3,
 118                         },
 119                 },
 120                 {
 121                         []float64{
 122                                 1, 2,
 123                                 3, 4,
 124                                 5, 6,
 125                         },
 126                         3, 2,
 127                         1, 6,
 128                         9.539392014169458,
 129                         &Dense{
 130                                 mat: blas64.General{
 131                                         Rows: 3, Cols: 2,
 132                                         Stride: 2,
 133                                         Data:   []float64{1, 2, 3, 4, 5, 6},
 134                                 },
 135                                 capRows: 3, capCols: 2,
 136                         },
 137                 },
 138         } {
 139                 m := NewDense(test.rows, test.cols, test.a)
 140                 rows, cols := m.Dims()
 141                 if rows != test.rows {
 142                         t.Errorf("unexpected number of rows for test %d: got: %d want: %d", i, rows, test.rows)
 143                 }
 144                 if cols != test.cols {
 145                         t.Errorf("unexpected number of cols for test %d: got: %d want: %d", i, cols, test.cols)
 146                 }
 147                 if min := Min(m); min != test.min {
 148                         t.Errorf("unexpected min for test %d: got: %v want: %v", i, min, test.min)
 149                 }
 150                 if max := Max(m); max != test.max {
 151                         t.Errorf("unexpected max for test %d: got: %v want: %v", i, max, test.max)
 152                 }
 153                 if fro := Norm(m, 2); math.Abs(Norm(m, 2)-test.fro) > 1e-14 {
 154                         t.Errorf("unexpected Frobenius norm for test %d: got: %v want: %v", i, fro, test.fro)
 155                 }
 156                 if !reflect.DeepEqual(m, test.mat) {
 157                         t.Errorf("unexpected matrix for test %d", i)
 158                 }
 159                 if !Equal(m, test.mat) {
 160                         t.Errorf("matrix does not equal expected matrix for test %d", i)
 161                 }
 162         }
 163 }
 164
 165 func TestAtSet(t *testing.T) {
 166         for test, af := range [][][]float64{
 167                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, // even
 168                 {{1, 2}, {4, 5}, {7, 8}},          // wide
 169                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, //skinny
 170         } {
 171                 m := NewDense(flatten(af))
 172                 rows, cols := m.Dims()
 173                 for i := 0; i < rows; i++ {
 174                         for j := 0; j < cols; j++ {
 175                                 if m.At(i, j) != af[i][j] {
 176                                         t.Errorf("unexpected value for At(%d, %d) for test %d: got: %v want: %v",
 177                                                 i, j, test, m.At(i, j), af[i][j])
 178                                 }
 179
 180                                 v := float64(i * j)
 181                                 m.Set(i, j, v)
 182                                 if m.At(i, j) != v {
 183                                         t.Errorf("unexpected value for At(%d, %d) after Set(%[1]d, %d, %v) for test %d: got: %v want: %[3]v",
 184                                                 i, j, v, test, m.At(i, j))
 185                                 }
 186                         }
 187                 }
 188                 // Check access out of bounds fails
 189                 for _, row := range []int{-1, rows, rows + 1} {
 190                         panicked, message := panics(func() { m.At(row, 0) })
 191                         if !panicked || message != ErrRowAccess.Error() {
 192                                 t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
 193                         }
 194                 }
 195                 for _, col := range []int{-1, cols, cols + 1} {
 196                         panicked, message := panics(func() { m.At(0, col) })
 197                         if !panicked || message != ErrColAccess.Error() {
 198                                 t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
 199                         }
 200                 }
 201
 202                 // Check Set out of bounds
 203                 for _, row := range []int{-1, rows, rows + 1} {
 204                         panicked, message := panics(func() { m.Set(row, 0, 1.2) })
 205                         if !panicked || message != ErrRowAccess.Error() {
 206                                 t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
 207                         }
 208                 }
 209                 for _, col := range []int{-1, cols, cols + 1} {
 210                         panicked, message := panics(func() { m.Set(0, col, 1.2) })
 211                         if !panicked || message != ErrColAccess.Error() {
 212                                 t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
 213                         }
 214                 }
 215         }
 216 }
 217
 218 func TestSetRowColumn(t *testing.T) {
 219         for _, as := range [][][]float64{
 220                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 221                 {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}},
 222                 {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}},
 223         } {
 224                 for ri, row := range as {
 225                         a := NewDense(flatten(as))
 226                         m := &Dense{}
 227                         m.Clone(a)
 228                         a.SetRow(ri, make([]float64, a.mat.Cols))
 229                         m.Sub(m, a)
 230                         nt := Norm(m, 2)
 231                         nr := floats.Norm(row, 2)
 232                         if math.Abs(nt-nr) > 1e-14 {
 233                                 t.Errorf("Row %d norm mismatch, want: %g, got: %g", ri, nr, nt)
 234                         }
 235                 }
 236
 237                 for ci := range as[0] {
 238                         a := NewDense(flatten(as))
 239                         m := &Dense{}
 240                         m.Clone(a)
 241                         a.SetCol(ci, make([]float64, a.mat.Rows))
 242                         col := make([]float64, a.mat.Rows)
 243                         for j := range col {
 244                                 col[j] = float64(ci + 1 + j*a.mat.Cols)
 245                         }
 246                         m.Sub(m, a)
 247                         nt := Norm(m, 2)
 248                         nc := floats.Norm(col, 2)
 249                         if math.Abs(nt-nc) > 1e-14 {
 250                                 t.Errorf("Column %d norm mismatch, want: %g, got: %g", ci, nc, nt)
 251                         }
 252                 }
 253         }
 254 }
 255
 256 func TestRowColView(t *testing.T) {
 257         for _, test := range []struct {
 258                 mat [][]float64
 259         }{
 260                 {
 261                         mat: [][]float64{
 262                                 {1, 2, 3, 4, 5},
 263                                 {6, 7, 8, 9, 10},
 264                                 {11, 12, 13, 14, 15},
 265                                 {16, 17, 18, 19, 20},
 266                                 {21, 22, 23, 24, 25},
 267                         },
 268                 },
 269                 {
 270                         mat: [][]float64{
 271                                 {1, 2, 3, 4},
 272                                 {6, 7, 8, 9},
 273                                 {11, 12, 13, 14},
 274                                 {16, 17, 18, 19},
 275                                 {21, 22, 23, 24},
 276                         },
 277                 },
 278                 {
 279                         mat: [][]float64{
 280                                 {1, 2, 3, 4, 5},
 281                                 {6, 7, 8, 9, 10},
 282                                 {11, 12, 13, 14, 15},
 283                                 {16, 17, 18, 19, 20},
 284                         },
 285                 },
 286         } {
 287                 // This over cautious approach to building a matrix data
 288                 // slice is to ensure that changes to flatten in the future
 289                 // do not mask a regression to the issue identified in
 290                 // gonum/matrix#110.
 291                 rows, cols, flat := flatten(test.mat)
 292                 m := NewDense(rows, cols, flat[:len(flat):len(flat)])
 293
 294                 for _, row := range []int{-1, rows, rows + 1} {
 295                         panicked, message := panics(func() { m.At(row, 0) })
 296                         if !panicked || message != ErrRowAccess.Error() {
 297                                 t.Errorf("expected panic for invalid row access rows=%d r=%d", rows, row)
 298                         }
 299                 }
 300                 for _, col := range []int{-1, cols, cols + 1} {
 301                         panicked, message := panics(func() { m.At(0, col) })
 302                         if !panicked || message != ErrColAccess.Error() {
 303                                 t.Errorf("expected panic for invalid column access cols=%d c=%d", cols, col)
 304                         }
 305                 }
 306
 307                 for i := 0; i < rows; i++ {
 308                         vr := m.RowView(i)
 309                         if vr.Len() != cols {
 310                                 t.Errorf("unexpected number of columns: got: %d want: %d", vr.Len(), cols)
 311                         }
 312                         for j := 0; j < cols; j++ {
 313                                 if got := vr.At(j, 0); got != test.mat[i][j] {
 314                                         t.Errorf("unexpected value for row.At(%d, 0): got: %v want: %v",
 315                                                 j, got, test.mat[i][j])
 316                                 }
 317                         }
 318                 }
 319                 for j := 0; j < cols; j++ {
 320                         vc := m.ColView(j)
 321                         if vc.Len() != rows {
 322                                 t.Errorf("unexpected number of rows: got: %d want: %d", vc.Len(), rows)
 323                         }
 324                         for i := 0; i < rows; i++ {
 325                                 if got := vc.At(i, 0); got != test.mat[i][j] {
 326                                         t.Errorf("unexpected value for col.At(%d, 0): got: %v want: %v",
 327                                                 i, got, test.mat[i][j])
 328                                 }
 329                         }
 330                 }
 331                 m = m.Slice(1, rows-1, 1, cols-1).(*Dense)
 332                 for i := 1; i < rows-1; i++ {
 333                         vr := m.RowView(i - 1)
 334                         if vr.Len() != cols-2 {
 335                                 t.Errorf("unexpected number of columns: got: %d want: %d", vr.Len(), cols-2)
 336                         }
 337                         for j := 1; j < cols-1; j++ {
 338                                 if got := vr.At(j-1, 0); got != test.mat[i][j] {
 339                                         t.Errorf("unexpected value for row.At(%d, 0): got: %v want: %v",
 340                                                 j-1, got, test.mat[i][j])
 341                                 }
 342                         }
 343                 }
 344                 for j := 1; j < cols-1; j++ {
 345                         vc := m.ColView(j - 1)
 346                         if vc.Len() != rows-2 {
 347                                 t.Errorf("unexpected number of rows: got: %d want: %d", vc.Len(), rows-2)
 348                         }
 349                         for i := 1; i < rows-1; i++ {
 350                                 if got := vc.At(i-1, 0); got != test.mat[i][j] {
 351                                         t.Errorf("unexpected value for col.At(%d, 0): got: %v want: %v",
 352                                                 i-1, got, test.mat[i][j])
 353                                 }
 354                         }
 355                 }
 356         }
 357 }
 358
 359 func TestGrow(t *testing.T) {
 360         m := &Dense{}
 361         m = m.Grow(10, 10).(*Dense)
 362         rows, cols := m.Dims()
 363         capRows, capCols := m.Caps()
 364         if rows != 10 {
 365                 t.Errorf("unexpected value for rows: got: %d want: 10", rows)
 366         }
 367         if cols != 10 {
 368                 t.Errorf("unexpected value for cols: got: %d want: 10", cols)
 369         }
 370         if capRows != 10 {
 371                 t.Errorf("unexpected value for capRows: got: %d want: 10", capRows)
 372         }
 373         if capCols != 10 {
 374                 t.Errorf("unexpected value for capCols: got: %d want: 10", capCols)
 375         }
 376
 377         // Test grow within caps is in-place.
 378         m.Set(1, 1, 1)
 379         v := m.Slice(1, 5, 1, 5).(*Dense)
 380         if v.At(0, 0) != m.At(1, 1) {
 381                 t.Errorf("unexpected viewed element value: got: %v want: %v", v.At(0, 0), m.At(1, 1))
 382         }
 383         v = v.Grow(5, 5).(*Dense)
 384         if !Equal(v, m.Slice(1, 10, 1, 10)) {
 385                 t.Error("unexpected view value after grow")
 386         }
 387
 388         // Test grow bigger than caps copies.
 389         v = v.Grow(5, 5).(*Dense)
 390         if !Equal(v.Slice(0, 9, 0, 9), m.Slice(1, 10, 1, 10)) {
 391                 t.Error("unexpected mismatched common view value after grow")
 392         }
 393         v.Set(0, 0, 0)
 394         if Equal(v.Slice(0, 9, 0, 9), m.Slice(1, 10, 1, 10)) {
 395                 t.Error("unexpected matching view value after grow past capacity")
 396         }
 397
 398         // Test grow uses existing data slice when matrix is zero size.
 399         v.Reset()
 400         p, l := &v.mat.Data[:1][0], cap(v.mat.Data)
 401         *p = 1 // This element is at position (-1, -1) relative to v and so should not be visible.
 402         v = v.Grow(5, 5).(*Dense)
 403         if &v.mat.Data[:1][0] != p {
 404                 t.Error("grow unexpectedly copied slice within cap limit")
 405         }
 406         if cap(v.mat.Data) != l {
 407                 t.Errorf("unexpected change in data slice capacity: got: %d want: %d", cap(v.mat.Data), l)
 408         }
 409         if v.At(0, 0) != 0 {
 410                 t.Errorf("unexpected value for At(0, 0): got: %v want: 0", v.At(0, 0))
 411         }
 412 }
 413
 414 func TestAdd(t *testing.T) {
 415         for i, test := range []struct {
 416                 a, b, r [][]float64
 417         }{
 418                 {
 419                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 420                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 421                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 422                 },
 423                 {
 424                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 425                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 426                         [][]float64{{2, 2, 2}, {2, 2, 2}, {2, 2, 2}},
 427                 },
 428                 {
 429                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 430                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 431                         [][]float64{{2, 0, 0}, {0, 2, 0}, {0, 0, 2}},
 432                 },
 433                 {
 434                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 435                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 436                         [][]float64{{-2, 0, 0}, {0, -2, 0}, {0, 0, -2}},
 437                 },
 438                 {
 439                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 440                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 441                         [][]float64{{2, 4, 6}, {8, 10, 12}},
 442                 },
 443         } {
 444                 a := NewDense(flatten(test.a))
 445                 b := NewDense(flatten(test.b))
 446                 r := NewDense(flatten(test.r))
 447
 448                 var temp Dense
 449                 temp.Add(a, b)
 450                 if !Equal(&temp, r) {
 451                         t.Errorf("unexpected result from Add for test %d %v Add %v: got: %v want: %v",
 452                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 453                 }
 454
 455                 zero(temp.mat.Data)
 456                 temp.Add(a, b)
 457                 if !Equal(&temp, r) {
 458                         t.Errorf("unexpected result from Add for test %d %v Add %v: got: %v want: %v",
 459                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 460                 }
 461
 462                 // These probably warrant a better check and failure. They should never happen in the wild though.
 463                 temp.mat.Data = nil
 464                 panicked, message := panics(func() { temp.Add(a, b) })
 465                 if !panicked || message != "runtime error: index out of range" {
 466                         t.Error("exected runtime panic for nil data slice")
 467                 }
 468
 469                 a.Add(a, b)
 470                 if !Equal(a, r) {
 471                         t.Errorf("unexpected result from Add for test %d %v Add %v: got: %v want: %v",
 472                                 i, test.a, test.b, unflatten(a.mat.Rows, a.mat.Cols, a.mat.Data), test.r)
 473                 }
 474         }
 475
 476         panicked, message := panics(func() {
 477                 m := NewDense(10, 10, nil)
 478                 a := NewDense(5, 5, nil)
 479                 m.Slice(1, 6, 1, 6).(*Dense).Add(a, m.Slice(2, 7, 2, 7))
 480         })
 481         if !panicked {
 482                 t.Error("expected panic for overlapping matrices")
 483         }
 484         if message != regionOverlap {
 485                 t.Errorf("unexpected panic message: got: %q want: %q", message, regionOverlap)
 486         }
 487
 488         method := func(receiver, a, b Matrix) {
 489                 type Adder interface {
 490                         Add(a, b Matrix)
 491                 }
 492                 rd := receiver.(Adder)
 493                 rd.Add(a, b)
 494         }
 495         denseComparison := func(receiver, a, b *Dense) {
 496                 receiver.Add(a, b)
 497         }
 498         testTwoInput(t, "Add", &Dense{}, method, denseComparison, legalTypesAll, legalSizeSameRectangular, 1e-14)
 499 }
 500
 501 func TestSub(t *testing.T) {
 502         for i, test := range []struct {
 503                 a, b, r [][]float64
 504         }{
 505                 {
 506                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 507                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 508                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 509                 },
 510                 {
 511                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 512                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 513                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 514                 },
 515                 {
 516                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 517                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 518                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 519                 },
 520                 {
 521                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 522                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 523                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 524                 },
 525                 {
 526                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 527                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 528                         [][]float64{{0, 0, 0}, {0, 0, 0}},
 529                 },
 530         } {
 531                 a := NewDense(flatten(test.a))
 532                 b := NewDense(flatten(test.b))
 533                 r := NewDense(flatten(test.r))
 534
 535                 var temp Dense
 536                 temp.Sub(a, b)
 537                 if !Equal(&temp, r) {
 538                         t.Errorf("unexpected result from Sub for test %d %v Sub %v: got: %v want: %v",
 539                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 540                 }
 541
 542                 zero(temp.mat.Data)
 543                 temp.Sub(a, b)
 544                 if !Equal(&temp, r) {
 545                         t.Errorf("unexpected result from Sub for test %d %v Sub %v: got: %v want: %v",
 546                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 547                 }
 548
 549                 // These probably warrant a better check and failure. They should never happen in the wild though.
 550                 temp.mat.Data = nil
 551                 panicked, message := panics(func() { temp.Sub(a, b) })
 552                 if !panicked || message != "runtime error: index out of range" {
 553                         t.Error("exected runtime panic for nil data slice")
 554                 }
 555
 556                 a.Sub(a, b)
 557                 if !Equal(a, r) {
 558                         t.Errorf("unexpected result from Sub for test %d %v Sub %v: got: %v want: %v",
 559                                 i, test.a, test.b, unflatten(a.mat.Rows, a.mat.Cols, a.mat.Data), test.r)
 560                 }
 561         }
 562
 563         panicked, message := panics(func() {
 564                 m := NewDense(10, 10, nil)
 565                 a := NewDense(5, 5, nil)
 566                 m.Slice(1, 6, 1, 6).(*Dense).Sub(a, m.Slice(2, 7, 2, 7))
 567         })
 568         if !panicked {
 569                 t.Error("expected panic for overlapping matrices")
 570         }
 571         if message != regionOverlap {
 572                 t.Errorf("unexpected panic message: got: %q want: %q", message, regionOverlap)
 573         }
 574
 575         method := func(receiver, a, b Matrix) {
 576                 type Suber interface {
 577                         Sub(a, b Matrix)
 578                 }
 579                 rd := receiver.(Suber)
 580                 rd.Sub(a, b)
 581         }
 582         denseComparison := func(receiver, a, b *Dense) {
 583                 receiver.Sub(a, b)
 584         }
 585         testTwoInput(t, "Sub", &Dense{}, method, denseComparison, legalTypesAll, legalSizeSameRectangular, 1e-14)
 586 }
 587
 588 func TestMulElem(t *testing.T) {
 589         for i, test := range []struct {
 590                 a, b, r [][]float64
 591         }{
 592                 {
 593                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 594                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 595                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 596                 },
 597                 {
 598                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 599                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 600                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 601                 },
 602                 {
 603                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 604                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 605                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 606                 },
 607                 {
 608                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 609                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 610                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 611                 },
 612                 {
 613                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 614                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 615                         [][]float64{{1, 4, 9}, {16, 25, 36}},
 616                 },
 617         } {
 618                 a := NewDense(flatten(test.a))
 619                 b := NewDense(flatten(test.b))
 620                 r := NewDense(flatten(test.r))
 621
 622                 var temp Dense
 623                 temp.MulElem(a, b)
 624                 if !Equal(&temp, r) {
 625                         t.Errorf("unexpected result from MulElem for test %d %v MulElem %v: got: %v want: %v",
 626                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 627                 }
 628
 629                 zero(temp.mat.Data)
 630                 temp.MulElem(a, b)
 631                 if !Equal(&temp, r) {
 632                         t.Errorf("unexpected result from MulElem for test %d %v MulElem %v: got: %v want: %v",
 633                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 634                 }
 635
 636                 // These probably warrant a better check and failure. They should never happen in the wild though.
 637                 temp.mat.Data = nil
 638                 panicked, message := panics(func() { temp.MulElem(a, b) })
 639                 if !panicked || message != "runtime error: index out of range" {
 640                         t.Error("exected runtime panic for nil data slice")
 641                 }
 642
 643                 a.MulElem(a, b)
 644                 if !Equal(a, r) {
 645                         t.Errorf("unexpected result from MulElem for test %d %v MulElem %v: got: %v want: %v",
 646                                 i, test.a, test.b, unflatten(a.mat.Rows, a.mat.Cols, a.mat.Data), test.r)
 647                 }
 648         }
 649
 650         panicked, message := panics(func() {
 651                 m := NewDense(10, 10, nil)
 652                 a := NewDense(5, 5, nil)
 653                 m.Slice(1, 6, 1, 6).(*Dense).MulElem(a, m.Slice(2, 7, 2, 7))
 654         })
 655         if !panicked {
 656                 t.Error("expected panic for overlapping matrices")
 657         }
 658         if message != regionOverlap {
 659                 t.Errorf("unexpected panic message: got: %q want: %q", message, regionOverlap)
 660         }
 661
 662         method := func(receiver, a, b Matrix) {
 663                 type ElemMuler interface {
 664                         MulElem(a, b Matrix)
 665                 }
 666                 rd := receiver.(ElemMuler)
 667                 rd.MulElem(a, b)
 668         }
 669         denseComparison := func(receiver, a, b *Dense) {
 670                 receiver.MulElem(a, b)
 671         }
 672         testTwoInput(t, "MulElem", &Dense{}, method, denseComparison, legalTypesAll, legalSizeSameRectangular, 1e-14)
 673 }
 674
 675 // A comparison that treats NaNs as equal, for testing.
 676 func (m *Dense) same(b Matrix) bool {
 677         br, bc := b.Dims()
 678         if br != m.mat.Rows || bc != m.mat.Cols {
 679                 return false
 680         }
 681         for r := 0; r < br; r++ {
 682                 for c := 0; c < bc; c++ {
 683                         if av, bv := m.At(r, c), b.At(r, c); av != bv && !(math.IsNaN(av) && math.IsNaN(bv)) {
 684                                 return false
 685                         }
 686                 }
 687         }
 688         return true
 689 }
 690
 691 func TestDivElem(t *testing.T) {
 692         for i, test := range []struct {
 693                 a, b, r [][]float64
 694         }{
 695                 {
 696                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 697                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 698                         [][]float64{{math.Inf(1), math.NaN(), math.NaN()}, {math.NaN(), math.Inf(1), math.NaN()}, {math.NaN(), math.NaN(), math.Inf(1)}},
 699                 },
 700                 {
 701                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 702                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 703                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 704                 },
 705                 {
 706                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 707                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 708                         [][]float64{{1, math.NaN(), math.NaN()}, {math.NaN(), 1, math.NaN()}, {math.NaN(), math.NaN(), 1}},
 709                 },
 710                 {
 711                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 712                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 713                         [][]float64{{1, math.NaN(), math.NaN()}, {math.NaN(), 1, math.NaN()}, {math.NaN(), math.NaN(), 1}},
 714                 },
 715                 {
 716                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 717                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 718                         [][]float64{{1, 1, 1}, {1, 1, 1}},
 719                 },
 720         } {
 721                 a := NewDense(flatten(test.a))
 722                 b := NewDense(flatten(test.b))
 723                 r := NewDense(flatten(test.r))
 724
 725                 var temp Dense
 726                 temp.DivElem(a, b)
 727                 if !temp.same(r) {
 728                         t.Errorf("unexpected result from DivElem for test %d %v DivElem %v: got: %v want: %v",
 729                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 730                 }
 731
 732                 zero(temp.mat.Data)
 733                 temp.DivElem(a, b)
 734                 if !temp.same(r) {
 735                         t.Errorf("unexpected result from DivElem for test %d %v DivElem %v: got: %v want: %v",
 736                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 737                 }
 738
 739                 // These probably warrant a better check and failure. They should never happen in the wild though.
 740                 temp.mat.Data = nil
 741                 panicked, message := panics(func() { temp.DivElem(a, b) })
 742                 if !panicked || message != "runtime error: index out of range" {
 743                         t.Error("exected runtime panic for nil data slice")
 744                 }
 745
 746                 a.DivElem(a, b)
 747                 if !a.same(r) {
 748                         t.Errorf("unexpected result from DivElem for test %d %v DivElem %v: got: %v want: %v",
 749                                 i, test.a, test.b, unflatten(a.mat.Rows, a.mat.Cols, a.mat.Data), test.r)
 750                 }
 751         }
 752
 753         panicked, message := panics(func() {
 754                 m := NewDense(10, 10, nil)
 755                 a := NewDense(5, 5, nil)
 756                 m.Slice(1, 6, 1, 6).(*Dense).DivElem(a, m.Slice(2, 7, 2, 7))
 757         })
 758         if !panicked {
 759                 t.Error("expected panic for overlapping matrices")
 760         }
 761         if message != regionOverlap {
 762                 t.Errorf("unexpected panic message: got: %q want: %q", message, regionOverlap)
 763         }
 764
 765         method := func(receiver, a, b Matrix) {
 766                 type ElemDiver interface {
 767                         DivElem(a, b Matrix)
 768                 }
 769                 rd := receiver.(ElemDiver)
 770                 rd.DivElem(a, b)
 771         }
 772         denseComparison := func(receiver, a, b *Dense) {
 773                 receiver.DivElem(a, b)
 774         }
 775         testTwoInput(t, "DivElem", &Dense{}, method, denseComparison, legalTypesAll, legalSizeSameRectangular, 1e-14)
 776 }
 777
 778 func TestMul(t *testing.T) {
 779         for i, test := range []struct {
 780                 a, b, r [][]float64
 781         }{
 782                 {
 783                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 784                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 785                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
 786                 },
 787                 {
 788                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 789                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
 790                         [][]float64{{3, 3, 3}, {3, 3, 3}, {3, 3, 3}},
 791                 },
 792                 {
 793                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 794                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 795                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 796                 },
 797                 {
 798                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 799                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
 800                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 801                 },
 802                 {
 803                         [][]float64{{1, 2, 3}, {4, 5, 6}},
 804                         [][]float64{{1, 2}, {3, 4}, {5, 6}},
 805                         [][]float64{{22, 28}, {49, 64}},
 806                 },
 807                 {
 808                         [][]float64{{0, 1, 1}, {0, 1, 1}, {0, 1, 1}},
 809                         [][]float64{{0, 1, 1}, {0, 1, 1}, {0, 1, 1}},
 810                         [][]float64{{0, 2, 2}, {0, 2, 2}, {0, 2, 2}},
 811                 },
 812         } {
 813                 a := NewDense(flatten(test.a))
 814                 b := NewDense(flatten(test.b))
 815                 r := NewDense(flatten(test.r))
 816
 817                 var temp Dense
 818                 temp.Mul(a, b)
 819                 if !Equal(&temp, r) {
 820                         t.Errorf("unexpected result from Mul for test %d %v Mul %v: got: %v want: %v",
 821                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 822                 }
 823
 824                 zero(temp.mat.Data)
 825                 temp.Mul(a, b)
 826                 if !Equal(&temp, r) {
 827                         t.Errorf("unexpected result from Mul for test %d %v Mul %v: got: %v want: %v",
 828                                 i, test.a, test.b, unflatten(temp.mat.Rows, temp.mat.Cols, temp.mat.Data), test.r)
 829                 }
 830
 831                 // These probably warrant a better check and failure. They should never happen in the wild though.
 832                 temp.mat.Data = nil
 833                 panicked, message := panics(func() { temp.Mul(a, b) })
 834                 if !panicked || message != "blas: index of c out of range" {
 835                         if message != "" {
 836                                 t.Errorf("expected runtime panic for nil data slice: got %q", message)
 837                         } else {
 838                                 t.Error("expected runtime panic for nil data slice")
 839                         }
 840                 }
 841         }
 842
 843         panicked, message := panics(func() {
 844                 m := NewDense(10, 10, nil)
 845                 a := NewDense(5, 5, nil)
 846                 m.Slice(1, 6, 1, 6).(*Dense).Mul(a, m.Slice(2, 7, 2, 7))
 847         })
 848         if !panicked {
 849                 t.Error("expected panic for overlapping matrices")
 850         }
 851         if message != regionOverlap {
 852                 t.Errorf("unexpected panic message: got: %q want: %q", message, regionOverlap)
 853         }
 854
 855         method := func(receiver, a, b Matrix) {
 856                 type Muler interface {
 857                         Mul(a, b Matrix)
 858                 }
 859                 rd := receiver.(Muler)
 860                 rd.Mul(a, b)
 861         }
 862         denseComparison := func(receiver, a, b *Dense) {
 863                 receiver.Mul(a, b)
 864         }
 865         legalSizeMul := func(ar, ac, br, bc int) bool {
 866                 return ac == br
 867         }
 868         testTwoInput(t, "Mul", &Dense{}, method, denseComparison, legalTypesAll, legalSizeMul, 1e-14)
 869 }
 870
 871 func randDense(size int, rho float64, rnd func() float64) (*Dense, error) {
 872         if size == 0 {
 873                 return nil, ErrZeroLength
 874         }
 875         d := &Dense{
 876                 mat: blas64.General{
 877                         Rows: size, Cols: size, Stride: size,
 878                         Data: make([]float64, size*size),
 879                 },
 880                 capRows: size, capCols: size,
 881         }
 882         for i := 0; i < size; i++ {
 883                 for j := 0; j < size; j++ {
 884                         if rand.Float64() < rho {
 885                                 d.Set(i, j, rnd())
 886                         }
 887                 }
 888         }
 889         return d, nil
 890 }
 891
 892 func TestExp(t *testing.T) {
 893         for i, test := range []struct {
 894                 a    [][]float64
 895                 want [][]float64
 896                 mod  func(*Dense)
 897         }{
 898                 {
 899                         a:    [][]float64{{-49, 24}, {-64, 31}},
 900                         want: [][]float64{{-0.7357587581474017, 0.5518190996594223}, {-1.4715175990917921, 1.103638240717339}},
 901                 },
 902                 {
 903                         a:    [][]float64{{-49, 24}, {-64, 31}},
 904                         want: [][]float64{{-0.7357587581474017, 0.5518190996594223}, {-1.4715175990917921, 1.103638240717339}},
 905                         mod: func(a *Dense) {
 906                                 d := make([]float64, 100)
 907                                 for i := range d {
 908                                         d[i] = math.NaN()
 909                                 }
 910                                 *a = *NewDense(10, 10, d).Slice(1, 3, 1, 3).(*Dense)
 911                         },
 912                 },
 913                 {
 914                         a:    [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 915                         want: [][]float64{{2.71828182845905, 0, 0}, {0, 2.71828182845905, 0}, {0, 0, 2.71828182845905}},
 916                 },
 917         } {
 918                 var got Dense
 919                 if test.mod != nil {
 920                         test.mod(&got)
 921                 }
 922                 got.Exp(NewDense(flatten(test.a)))
 923                 if !EqualApprox(&got, NewDense(flatten(test.want)), 1e-12) {
 924                         t.Errorf("unexpected result for Exp test %d", i)
 925                 }
 926         }
 927 }
 928
 929 func TestPow(t *testing.T) {
 930         for i, test := range []struct {
 931                 a    [][]float64
 932                 n    int
 933                 mod  func(*Dense)
 934                 want [][]float64
 935         }{
 936                 {
 937                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 938                         n:    0,
 939                         want: [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 940                 },
 941                 {
 942                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 943                         n:    0,
 944                         want: [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
 945                         mod: func(a *Dense) {
 946                                 d := make([]float64, 100)
 947                                 for i := range d {
 948                                         d[i] = math.NaN()
 949                                 }
 950                                 *a = *NewDense(10, 10, d).Slice(1, 4, 1, 4).(*Dense)
 951                         },
 952                 },
 953                 {
 954                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 955                         n:    1,
 956                         want: [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 957                 },
 958                 {
 959                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 960                         n:    1,
 961                         want: [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 962                         mod: func(a *Dense) {
 963                                 d := make([]float64, 100)
 964                                 for i := range d {
 965                                         d[i] = math.NaN()
 966                                 }
 967                                 *a = *NewDense(10, 10, d).Slice(1, 4, 1, 4).(*Dense)
 968                         },
 969                 },
 970                 {
 971                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 972                         n:    2,
 973                         want: [][]float64{{30, 36, 42}, {66, 81, 96}, {102, 126, 150}},
 974                 },
 975                 {
 976                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 977                         n:    2,
 978                         want: [][]float64{{30, 36, 42}, {66, 81, 96}, {102, 126, 150}},
 979                         mod: func(a *Dense) {
 980                                 d := make([]float64, 100)
 981                                 for i := range d {
 982                                         d[i] = math.NaN()
 983                                 }
 984                                 *a = *NewDense(10, 10, d).Slice(1, 4, 1, 4).(*Dense)
 985                         },
 986                 },
 987                 {
 988                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 989                         n:    3,
 990                         want: [][]float64{{468, 576, 684}, {1062, 1305, 1548}, {1656, 2034, 2412}},
 991                 },
 992                 {
 993                         a:    [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 994                         n:    3,
 995                         want: [][]float64{{468, 576, 684}, {1062, 1305, 1548}, {1656, 2034, 2412}},
 996                         mod: func(a *Dense) {
 997                                 d := make([]float64, 100)
 998                                 for i := range d {
 999                                         d[i] = math.NaN()
1000                                 }
1001                                 *a = *NewDense(10, 10, d).Slice(1, 4, 1, 4).(*Dense)
1002                         },
1003                 },
1004         } {
1005                 var got Dense
1006                 if test.mod != nil {
1007                         test.mod(&got)
1008                 }
1009                 got.Pow(NewDense(flatten(test.a)), test.n)
1010                 if !EqualApprox(&got, NewDense(flatten(test.want)), 1e-12) {
1011                         t.Errorf("unexpected result for Pow test %d", i)
1012                 }
1013         }
1014 }
1015
1016 func TestScale(t *testing.T) {
1017         for _, f := range []float64{0.5, 1, 3} {
1018                 method := func(receiver, a Matrix) {
1019                         type Scaler interface {
1020                                 Scale(f float64, a Matrix)
1021                         }
1022                         rd := receiver.(Scaler)
1023                         rd.Scale(f, a)
1024                 }
1025                 denseComparison := func(receiver, a *Dense) {
1026                         receiver.Scale(f, a)
1027                 }
1028                 testOneInput(t, "Scale", &Dense{}, method, denseComparison, isAnyType, isAnySize, 1e-14)
1029         }
1030 }
1031
1032 func TestPowN(t *testing.T) {
1033         for i, test := range []struct {
1034                 a   [][]float64
1035                 mod func(*Dense)
1036         }{
1037                 {
1038                         a: [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
1039                 },
1040                 {
1041                         a: [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
1042                         mod: func(a *Dense) {
1043                                 d := make([]float64, 100)
1044                                 for i := range d {
1045                                         d[i] = math.NaN()
1046                                 }
1047                                 *a = *NewDense(10, 10, d).Slice(1, 4, 1, 4).(*Dense)
1048                         },
1049                 },
1050         } {
1051                 for n := 1; n <= 14; n++ {
1052                         var got, want Dense
1053                         if test.mod != nil {
1054                                 test.mod(&got)
1055                         }
1056                         got.Pow(NewDense(flatten(test.a)), n)
1057                         want.iterativePow(NewDense(flatten(test.a)), n)
1058                         if !Equal(&got, &want) {
1059                                 t.Errorf("unexpected result for iterative Pow test %d", i)
1060                         }
1061                 }
1062         }
1063 }
1064
1065 func (m *Dense) iterativePow(a Matrix, n int) {
1066         m.Clone(a)
1067         for i := 1; i < n; i++ {
1068                 m.Mul(m, a)
1069         }
1070 }
1071
1072 func TestCloneT(t *testing.T) {
1073         for i, test := range []struct {
1074                 a, want [][]float64
1075         }{
1076                 {
1077                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1078                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1079                 },
1080                 {
1081                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1082                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1083                 },
1084                 {
1085                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1086                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1087                 },
1088                 {
1089                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
1090                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
1091                 },
1092                 {
1093                         [][]float64{{1, 2, 3}, {4, 5, 6}},
1094                         [][]float64{{1, 4}, {2, 5}, {3, 6}},
1095                 },
1096         } {
1097                 a := NewDense(flatten(test.a))
1098                 want := NewDense(flatten(test.want))
1099
1100                 var got, gotT Dense
1101
1102                 for j := 0; j < 2; j++ {
1103                         got.Clone(a.T())
1104                         if !Equal(&got, want) {
1105                                 t.Errorf("expected transpose for test %d iteration %d: %v transpose = %v",
1106                                         i, j, test.a, test.want)
1107                         }
1108                         gotT.Clone(got.T())
1109                         if !Equal(&gotT, a) {
1110                                 t.Errorf("expected transpose for test %d iteration %d: %v transpose = %v",
1111                                         i, j, test.a, test.want)
1112                         }
1113
1114                         zero(got.mat.Data)
1115                 }
1116         }
1117 }
1118
1119 func TestCopyT(t *testing.T) {
1120         for i, test := range []struct {
1121                 a, want [][]float64
1122         }{
1123                 {
1124                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1125                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1126                 },
1127                 {
1128                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1129                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1130                 },
1131                 {
1132                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1133                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1134                 },
1135                 {
1136                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
1137                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
1138                 },
1139                 {
1140                         [][]float64{{1, 2, 3}, {4, 5, 6}},
1141                         [][]float64{{1, 4}, {2, 5}, {3, 6}},
1142                 },
1143         } {
1144                 a := NewDense(flatten(test.a))
1145                 want := NewDense(flatten(test.want))
1146
1147                 ar, ac := a.Dims()
1148                 got := NewDense(ac, ar, nil)
1149                 rr := NewDense(ar, ac, nil)
1150
1151                 for j := 0; j < 2; j++ {
1152                         got.Copy(a.T())
1153                         if !Equal(got, want) {
1154                                 t.Errorf("expected transpose for test %d iteration %d: %v transpose = %v",
1155                                         i, j, test.a, test.want)
1156                         }
1157                         rr.Copy(got.T())
1158                         if !Equal(rr, a) {
1159                                 t.Errorf("expected transpose for test %d iteration %d: %v transpose = %v",
1160                                         i, j, test.a, test.want)
1161                         }
1162
1163                         zero(got.mat.Data)
1164                 }
1165         }
1166 }
1167
1168 func TestCopyDenseAlias(t *testing.T) {
1169         for _, trans := range []bool{false, true} {
1170                 for di := 0; di < 2; di++ {
1171                         for dj := 0; dj < 2; dj++ {
1172                                 for si := 0; si < 2; si++ {
1173                                         for sj := 0; sj < 2; sj++ {
1174                                                 a := NewDense(3, 3, []float64{
1175                                                         1, 2, 3,
1176                                                         4, 5, 6,
1177                                                         7, 8, 9,
1178                                                 })
1179                                                 src := a.Slice(si, si+2, sj, sj+2)
1180                                                 want := DenseCopyOf(src)
1181                                                 got := a.Slice(di, di+2, dj, dj+2).(*Dense)
1182
1183                                                 if trans {
1184                                                         panicked, _ := panics(func() { got.Copy(src.T()) })
1185                                                         if !panicked {
1186                                                                 t.Errorf("expected panic for transpose aliased copy with offsets dst(%d,%d) src(%d,%d):\ngot:\n%v\nwant:\n%v",
1187                                                                         di, dj, si, sj, Formatted(got), Formatted(want),
1188                                                                 )
1189                                                         }
1190                                                         continue
1191                                                 }
1192
1193                                                 got.Copy(src)
1194                                                 if !Equal(got, want) {
1195                                                         t.Errorf("unexpected aliased copy result with offsets dst(%d,%d) src(%d,%d):\ngot:\n%v\nwant:\n%v",
1196                                                                 di, dj, si, sj, Formatted(got), Formatted(want),
1197                                                         )
1198                                                 }
1199                                         }
1200                                 }
1201                         }
1202                 }
1203         }
1204 }
1205
1206 func TestCopyVecDenseAlias(t *testing.T) {
1207         for _, horiz := range []bool{false, true} {
1208                 for do := 0; do < 2; do++ {
1209                         for di := 0; di < 3; di++ {
1210                                 for si := 0; si < 3; si++ {
1211                                         a := NewDense(3, 3, []float64{
1212                                                 1, 2, 3,
1213                                                 4, 5, 6,
1214                                                 7, 8, 9,
1215                                         })
1216                                         var src Vector
1217                                         var want *Dense
1218                                         if horiz {
1219                                                 src = a.RowView(si)
1220                                                 want = DenseCopyOf(a.Slice(si, si+1, 0, 2))
1221                                         } else {
1222                                                 src = a.ColView(si)
1223                                                 want = DenseCopyOf(a.Slice(0, 2, si, si+1))
1224                                         }
1225
1226                                         var got *Dense
1227                                         if horiz {
1228                                                 got = a.Slice(di, di+1, do, do+2).(*Dense)
1229                                                 got.Copy(src.T())
1230                                         } else {
1231                                                 got = a.Slice(do, do+2, di, di+1).(*Dense)
1232                                                 got.Copy(src)
1233                                         }
1234
1235                                         if !Equal(got, want) {
1236                                                 t.Errorf("unexpected aliased copy result with offsets dst(%d) src(%d):\ngot:\n%v\nwant:\n%v",
1237                                                         di, si, Formatted(got), Formatted(want),
1238                                                 )
1239                                         }
1240                                 }
1241                         }
1242                 }
1243         }
1244 }
1245
1246 func identity(r, c int, v float64) float64 { return v }
1247
1248 func TestApply(t *testing.T) {
1249         for i, test := range []struct {
1250                 a, want [][]float64
1251                 fn      func(r, c int, v float64) float64
1252         }{
1253                 {
1254                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1255                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1256                         identity,
1257                 },
1258                 {
1259                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1260                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1261                         identity,
1262                 },
1263                 {
1264                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1265                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1266                         identity,
1267                 },
1268                 {
1269                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
1270                         [][]float64{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
1271                         identity,
1272                 },
1273                 {
1274                         [][]float64{{1, 2, 3}, {4, 5, 6}},
1275                         [][]float64{{1, 2, 3}, {4, 5, 6}},
1276                         identity,
1277                 },
1278                 {
1279                         [][]float64{{1, 2, 3}, {4, 5, 6}},
1280                         [][]float64{{2, 4, 6}, {8, 10, 12}},
1281                         func(r, c int, v float64) float64 { return v * 2 },
1282                 },
1283                 {
1284                         [][]float64{{1, 2, 3}, {4, 5, 6}},
1285                         [][]float64{{0, 2, 0}, {0, 5, 0}},
1286                         func(r, c int, v float64) float64 {
1287                                 if c == 1 {
1288                                         return v
1289                                 }
1290                                 return 0
1291                         },
1292                 },
1293                 {
1294                         [][]float64{{1, 2, 3}, {4, 5, 6}},
1295                         [][]float64{{0, 0, 0}, {4, 5, 6}},
1296                         func(r, c int, v float64) float64 {
1297                                 if r == 1 {
1298                                         return v
1299                                 }
1300                                 return 0
1301                         },
1302                 },
1303         } {
1304                 a := NewDense(flatten(test.a))
1305                 want := NewDense(flatten(test.want))
1306
1307                 var got Dense
1308
1309                 for j := 0; j < 2; j++ {
1310                         got.Apply(test.fn, a)
1311                         if !Equal(&got, want) {
1312                                 t.Errorf("unexpected result for test %d iteration %d: got: %v want: %v", i, j, got.mat.Data, want.mat.Data)
1313                         }
1314                 }
1315         }
1316
1317         for _, fn := range []func(r, c int, v float64) float64{
1318                 identity,
1319                 func(r, c int, v float64) float64 {
1320                         if r < c {
1321                                 return v
1322                         }
1323                         return -v
1324                 },
1325                 func(r, c int, v float64) float64 {
1326                         if r%2 == 0 && c%2 == 0 {
1327                                 return v
1328                         }
1329                         return -v
1330                 },
1331                 func(_, _ int, v float64) float64 { return v * v },
1332                 func(_, _ int, v float64) float64 { return -v },
1333         } {
1334                 method := func(receiver, x Matrix) {
1335                         type Applier interface {
1336                                 Apply(func(r, c int, v float64) float64, Matrix)
1337                         }
1338                         rd := receiver.(Applier)
1339                         rd.Apply(fn, x)
1340                 }
1341                 denseComparison := func(receiver, x *Dense) {
1342                         receiver.Apply(fn, x)
1343                 }
1344                 testOneInput(t, "Apply", &Dense{}, method, denseComparison, isAnyType, isAnySize, 0)
1345         }
1346 }
1347
1348 func TestClone(t *testing.T) {
1349         for i, test := range []struct {
1350                 a    [][]float64
1351                 i, j int
1352                 v    float64
1353         }{
1354                 {
1355                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1356                         1, 1,
1357                         1,
1358                 },
1359                 {
1360                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1361                         0, 0,
1362                         0,
1363                 },
1364         } {
1365                 a := NewDense(flatten(test.a))
1366                 b := *a
1367                 a.Clone(a)
1368                 a.Set(test.i, test.j, test.v)
1369
1370                 if Equal(&b, a) {
1371                         t.Errorf("unexpected mirror of write to cloned matrix for test %d: %v cloned and altered = %v",
1372                                 i, a, &b)
1373                 }
1374         }
1375 }
1376
1377 // TODO(kortschak) Roll this into testOneInput when it exists.
1378 func TestCopyPanic(t *testing.T) {
1379         for _, a := range []*Dense{
1380                 {},
1381                 {mat: blas64.General{Rows: 1}},
1382                 {mat: blas64.General{Cols: 1}},
1383         } {
1384                 var rows, cols int
1385                 m := NewDense(1, 1, nil)
1386                 panicked, message := panics(func() { rows, cols = m.Copy(a) })
1387                 if panicked {
1388                         t.Errorf("unexpected panic: %v", message)
1389                 }
1390                 if rows != 0 {
1391                         t.Errorf("unexpected rows: got: %d want: 0", rows)
1392                 }
1393                 if cols != 0 {
1394                         t.Errorf("unexpected cols: got: %d want: 0", cols)
1395                 }
1396         }
1397 }
1398
1399 func TestStack(t *testing.T) {
1400         for i, test := range []struct {
1401                 a, b, e [][]float64
1402         }{
1403                 {
1404                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1405                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1406                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1407                 },
1408                 {
1409                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1410                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1411                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}, {1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1412                 },
1413                 {
1414                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1415                         [][]float64{{0, 1, 0}, {0, 0, 1}, {1, 0, 0}},
1416                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 1, 0}, {0, 0, 1}, {1, 0, 0}},
1417                 },
1418         } {
1419                 a := NewDense(flatten(test.a))
1420                 b := NewDense(flatten(test.b))
1421
1422                 var s Dense
1423                 s.Stack(a, b)
1424
1425                 if !Equal(&s, NewDense(flatten(test.e))) {
1426                         t.Errorf("unexpected result for Stack test %d: %v stack %v = %v", i, a, b, s)
1427                 }
1428         }
1429
1430         method := func(receiver, a, b Matrix) {
1431                 type Stacker interface {
1432                         Stack(a, b Matrix)
1433                 }
1434                 rd := receiver.(Stacker)
1435                 rd.Stack(a, b)
1436         }
1437         denseComparison := func(receiver, a, b *Dense) {
1438                 receiver.Stack(a, b)
1439         }
1440         testTwoInput(t, "Stack", &Dense{}, method, denseComparison, legalTypesAll, legalSizeSameWidth, 0)
1441 }
1442
1443 func TestAugment(t *testing.T) {
1444         for i, test := range []struct {
1445                 a, b, e [][]float64
1446         }{
1447                 {
1448                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1449                         [][]float64{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1450                         [][]float64{{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}},
1451                 },
1452                 {
1453                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1454                         [][]float64{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}},
1455                         [][]float64{{1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1}},
1456                 },
1457                 {
1458                         [][]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
1459                         [][]float64{{0, 1, 0}, {0, 0, 1}, {1, 0, 0}},
1460                         [][]float64{{1, 0, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 1}, {0, 0, 1, 1, 0, 0}},
1461                 },
1462         } {
1463                 a := NewDense(flatten(test.a))
1464                 b := NewDense(flatten(test.b))
1465
1466                 var s Dense
1467                 s.Augment(a, b)
1468
1469                 if !Equal(&s, NewDense(flatten(test.e))) {
1470                         t.Errorf("unexpected result for Augment test %d: %v augment %v = %v", i, a, b, s)
1471                 }
1472         }
1473
1474         method := func(receiver, a, b Matrix) {
1475                 type Augmenter interface {
1476                         Augment(a, b Matrix)
1477                 }
1478                 rd := receiver.(Augmenter)
1479                 rd.Augment(a, b)
1480         }
1481         denseComparison := func(receiver, a, b *Dense) {
1482                 receiver.Augment(a, b)
1483         }
1484         testTwoInput(t, "Augment", &Dense{}, method, denseComparison, legalTypesAll, legalSizeSameHeight, 0)
1485 }
1486
1487 func TestRankOne(t *testing.T) {
1488         for i, test := range []struct {
1489                 x     []float64
1490                 y     []float64
1491                 m     [][]float64
1492                 alpha float64
1493         }{
1494                 {
1495                         x:     []float64{5},
1496                         y:     []float64{10},
1497                         m:     [][]float64{{2}},
1498                         alpha: -3,
1499                 },
1500                 {
1501                         x:     []float64{5, 6, 1},
1502                         y:     []float64{10},
1503                         m:     [][]float64{{2}, {-3}, {5}},
1504                         alpha: -3,
1505                 },
1506
1507                 {
1508                         x:     []float64{5},
1509                         y:     []float64{10, 15, 8},
1510                         m:     [][]float64{{2, -3, 5}},
1511                         alpha: -3,
1512                 },
1513                 {
1514                         x: []float64{1, 5},
1515                         y: []float64{10, 15},
1516                         m: [][]float64{
1517                                 {2, -3},
1518                                 {4, -1},
1519                         },
1520                         alpha: -3,
1521                 },
1522                 {
1523                         x: []float64{2, 3, 9},
1524                         y: []float64{8, 9},
1525                         m: [][]float64{
1526                                 {2, 3},
1527                                 {4, 5},
1528                                 {6, 7},
1529                         },
1530                         alpha: -3,
1531                 },
1532                 {
1533                         x: []float64{2, 3},
1534                         y: []float64{8, 9, 9},
1535                         m: [][]float64{
1536                                 {2, 3, 6},
1537                                 {4, 5, 7},
1538                         },
1539                         alpha: -3,
1540                 },
1541         } {
1542                 want := &Dense{}
1543                 xm := NewDense(len(test.x), 1, test.x)
1544                 ym := NewDense(1, len(test.y), test.y)
1545
1546                 want.Mul(xm, ym)
1547                 want.Scale(test.alpha, want)
1548                 want.Add(want, NewDense(flatten(test.m)))
1549
1550                 a := NewDense(flatten(test.m))
1551                 m := &Dense{}
1552                 // Check with a new matrix
1553                 m.RankOne(a, test.alpha, NewVecDense(len(test.x), test.x), NewVecDense(len(test.y), test.y))
1554                 if !Equal(m, want) {
1555                         t.Errorf("unexpected result for RankOne test %d iteration 0: got: %+v want: %+v", i, m, want)
1556                 }
1557                 // Check with the same matrix
1558                 a.RankOne(a, test.alpha, NewVecDense(len(test.x), test.x), NewVecDense(len(test.y), test.y))
1559                 if !Equal(a, want) {
1560                         t.Errorf("unexpected result for RankOne test %d iteration 1: got: %+v want: %+v", i, m, want)
1561                 }
1562         }
1563 }
1564
1565 func TestOuter(t *testing.T) {
1566         for i, test := range []struct {
1567                 x []float64
1568                 y []float64
1569         }{
1570                 {
1571                         x: []float64{5},
1572                         y: []float64{10},
1573                 },
1574                 {
1575                         x: []float64{5, 6, 1},
1576                         y: []float64{10},
1577                 },
1578
1579                 {
1580                         x: []float64{5},
1581                         y: []float64{10, 15, 8},
1582                 },
1583                 {
1584                         x: []float64{1, 5},
1585                         y: []float64{10, 15},
1586                 },
1587                 {
1588                         x: []float64{2, 3, 9},
1589                         y: []float64{8, 9},
1590                 },
1591                 {
1592                         x: []float64{2, 3},
1593                         y: []float64{8, 9, 9},
1594                 },
1595         } {
1596                 for _, f := range []float64{0.5, 1, 3} {
1597                         want := &Dense{}
1598                         xm := NewDense(len(test.x), 1, test.x)
1599                         ym := NewDense(1, len(test.y), test.y)
1600
1601                         want.Mul(xm, ym)
1602                         want.Scale(f, want)
1603
1604                         var m Dense
1605                         for j := 0; j < 2; j++ {
1606                                 // Check with a new matrix - and then again.
1607                                 m.Outer(f, NewVecDense(len(test.x), test.x), NewVecDense(len(test.y), test.y))
1608                                 if !Equal(&m, want) {
1609                                         t.Errorf("unexpected result for Outer test %d iteration %d scale %v: got: %+v want: %+v", i, j, f, m, want)
1610                                 }
1611                         }
1612                 }
1613         }
1614
1615         for _, alpha := range []float64{0, 1, -1, 2.3, -2.3} {
1616                 method := func(receiver, x, y Matrix) {
1617                         type outerer interface {
1618                                 Outer(alpha float64, x, y Vector)
1619                         }
1620                         m := receiver.(outerer)
1621                         m.Outer(alpha, x.(Vector), y.(Vector))
1622                 }
1623                 denseComparison := func(receiver, x, y *Dense) {
1624                         receiver.Mul(x, y.T())
1625                         receiver.Scale(alpha, receiver)
1626                 }
1627                 testTwoInput(t, "Outer", &Dense{}, method, denseComparison, legalTypesVectorVector, legalSizeVector, 1e-12)
1628         }
1629 }
1630
1631 func TestInverse(t *testing.T) {
1632         for i, test := range []struct {
1633                 a    Matrix
1634                 want Matrix // nil indicates that a is singular.
1635                 tol  float64
1636         }{
1637                 {
1638                         a: NewDense(3, 3, []float64{
1639                                 8, 1, 6,
1640                                 3, 5, 7,
1641                                 4, 9, 2,
1642                         }),
1643                         want: NewDense(3, 3, []float64{
1644                                 0.147222222222222, -0.144444444444444, 0.063888888888889,
1645                                 -0.061111111111111, 0.022222222222222, 0.105555555555556,
1646                                 -0.019444444444444, 0.188888888888889, -0.102777777777778,
1647                         }),
1648                         tol: 1e-14,
1649                 },
1650                 {
1651                         a: NewDense(3, 3, []float64{
1652                                 8, 1, 6,
1653                                 3, 5, 7,
1654                                 4, 9, 2,
1655                         }).T(),
1656                         want: NewDense(3, 3, []float64{
1657                                 0.147222222222222, -0.144444444444444, 0.063888888888889,
1658                                 -0.061111111111111, 0.022222222222222, 0.105555555555556,
1659                                 -0.019444444444444, 0.188888888888889, -0.102777777777778,
1660                         }).T(),
1661                         tol: 1e-14,
1662                 },
1663
1664                 // This case does not fail, but we do not guarantee that. The success
1665                 // is because the receiver and the input are aligned in the call to
1666                 // inverse. If there was a misalignment, the result would likely be
1667                 // incorrect and no shadowing panic would occur.
1668                 {
1669                         a: asBasicMatrix(NewDense(3, 3, []float64{
1670                                 8, 1, 6,
1671                                 3, 5, 7,
1672                                 4, 9, 2,
1673                         })),
1674                         want: NewDense(3, 3, []float64{
1675                                 0.147222222222222, -0.144444444444444, 0.063888888888889,
1676                                 -0.061111111111111, 0.022222222222222, 0.105555555555556,
1677                                 -0.019444444444444, 0.188888888888889, -0.102777777777778,
1678                         }),
1679                         tol: 1e-14,
1680                 },
1681
1682                 // The following case fails as it does not follow the shadowing rules.
1683                 // Specifically, the test extracts the underlying *Dense, and uses
1684                 // it as a receiver with the basicMatrix as input. The basicMatrix type
1685                 // allows shadowing of the input data without providing the Raw method
1686                 // required for detection of shadowing.
1687                 //
1688                 // We specifically state we do not check this case.
1689                 //
1690                 // {
1691                 //      a: asBasicMatrix(NewDense(3, 3, []float64{
1692                 //              8, 1, 6,
1693                 //              3, 5, 7,
1694                 //              4, 9, 2,
1695                 //      })).T(),
1696                 //      want: NewDense(3, 3, []float64{
1697                 //              0.147222222222222, -0.144444444444444, 0.063888888888889,
1698                 //              -0.061111111111111, 0.022222222222222, 0.105555555555556,
1699                 //              -0.019444444444444, 0.188888888888889, -0.102777777777778,
1700                 //      }).T(),
1701                 //      tol: 1e-14,
1702                 // },
1703
1704                 {
1705                         a: NewDense(4, 4, []float64{
1706                                 5, 2, 8, 7,
1707                                 4, 5, 8, 2,
1708                                 8, 5, 3, 2,
1709                                 8, 7, 7, 5,
1710                         }),
1711                         want: NewDense(4, 4, []float64{
1712                                 0.100548446069470, 0.021937842778793, 0.334552102376599, -0.283363802559415,
1713                                 -0.226691042047532, -0.067641681901280, -0.281535648994515, 0.457038391224863,
1714                                 0.080438756855576, 0.217550274223035, 0.067641681901280, -0.226691042047532,
1715                                 0.043875685557587, -0.244972577696527, -0.235831809872029, 0.330895795246801,
1716                         }),
1717                         tol: 1e-14,
1718                 },
1719
1720                 // Tests with singular matrix.
1721                 {
1722                         a: NewDense(1, 1, []float64{
1723                                 0,
1724                         }),
1725                 },
1726                 {
1727                         a: NewDense(2, 2, []float64{
1728                                 0, 0,
1729                                 0, 0,
1730                         }),
1731                 },
1732                 {
1733                         a: NewDense(2, 2, []float64{
1734                                 0, 0,
1735                                 0, 1,
1736                         }),
1737                 },
1738                 {
1739                         a: NewDense(3, 3, []float64{
1740                                 0, 0, 0,
1741                                 0, 0, 0,
1742                                 0, 0, 0,
1743                         }),
1744                 },
1745                 {
1746                         a: NewDense(4, 4, []float64{
1747                                 0, 0, 0, 0,
1748                                 0, 0, 0, 0,
1749                                 0, 0, 0, 0,
1750                                 0, 0, 0, 0,
1751                         }),
1752                 },
1753                 {
1754                         a: NewDense(4, 4, []float64{
1755                                 0, 0, 0, 0,
1756                                 0, 0, 0, 0,
1757                                 0, 0, 20, 20,
1758                                 0, 0, 20, 20,
1759                         }),
1760                 },
1761                 {
1762                         a: NewDense(4, 4, []float64{
1763                                 0, 1, 0, 0,
1764                                 0, 0, 1, 0,
1765                                 0, 0, 0, 1,
1766                                 0, 0, 0, 0,
1767                         }),
1768                 },
1769                 {
1770                         a: NewDense(4, 4, []float64{
1771                                 1, 1, 1, 1,
1772                                 1, 1, 1, 1,
1773                                 1, 1, 1, 1,
1774                                 1, 1, 1, 1,
1775                         }),
1776                 },
1777                 {
1778                         a: NewDense(5, 5, []float64{
1779                                 0, 1, 0, 0, 0,
1780                                 4, 0, 2, 0, 0,
1781                                 0, 3, 0, 3, 0,
1782                                 0, 0, 2, 0, 4,
1783                                 0, 0, 0, 1, 0,
1784                         }),
1785                 },
1786                 {
1787                         a: NewDense(5, 5, []float64{
1788                                 4, -1, -1, -1, -1,
1789                                 -1, 4, -1, -1, -1,
1790                                 -1, -1, 4, -1, -1,
1791                                 -1, -1, -1, 4, -1,
1792                                 -1, -1, -1, -1, 4,
1793                         }),
1794                 },
1795                 {
1796                         a: NewDense(5, 5, []float64{
1797                                 2, -1, 0, 0, -1,
1798                                 -1, 2, -1, 0, 0,
1799                                 0, -1, 2, -1, 0,
1800                                 0, 0, -1, 2, -1,
1801                                 -1, 0, 0, -1, 2,
1802                         }),
1803                 },
1804                 {
1805                         a: NewDense(5, 5, []float64{
1806                                 1, 2, 3, 5, 8,
1807                                 2, 3, 5, 8, 13,
1808                                 3, 5, 8, 13, 21,
1809                                 5, 8, 13, 21, 34,
1810                                 8, 13, 21, 34, 55,
1811                         }),
1812                 },
1813                 {
1814                         a: NewDense(8, 8, []float64{
1815                                 611, 196, -192, 407, -8, -52, -49, 29,
1816                                 196, 899, 113, -192, -71, -43, -8, -44,
1817                                 -192, 113, 899, 196, 61, 49, 8, 52,
1818                                 407, -192, 196, 611, 8, 44, 59, -23,
1819                                 -8, -71, 61, 8, 411, -599, 208, 208,
1820                                 -52, -43, 49, 44, -599, 411, 208, 208,
1821                                 -49, -8, 8, 59, 208, 208, 99, -911,
1822                                 29, -44, 52, -23, 208, 208, -911, 99,
1823                         }),
1824                 },
1825         } {
1826                 var got Dense
1827                 err := got.Inverse(test.a)
1828                 if test.want == nil {
1829                         if err == nil {
1830                                 t.Errorf("Case %d: expected error for singular matrix", i)
1831                         }
1832                         continue
1833                 }
1834                 if err != nil {
1835                         t.Errorf("Case %d: unexpected error: %v", i, err)
1836                         continue
1837                 }
1838                 if !equalApprox(&got, test.want, test.tol, false) {
1839                         t.Errorf("Case %d, inverse mismatch.", i)
1840                 }
1841                 var m Dense
1842                 m.Mul(&got, test.a)
1843                 r, _ := test.a.Dims()
1844                 d := make([]float64, r*r)
1845                 for i := 0; i < r*r; i += r + 1 {
1846                         d[i] = 1
1847                 }
1848                 eye := NewDense(r, r, d)
1849                 if !equalApprox(eye, &m, 1e-14, false) {
1850                         t.Errorf("Case %d, A^-1 * A != I", i)
1851                 }
1852
1853                 var tmp Dense
1854                 tmp.Clone(test.a)
1855                 aU, transposed := untranspose(test.a)
1856                 if transposed {
1857                         switch aU := aU.(type) {
1858                         case *Dense:
1859                                 err = aU.Inverse(test.a)
1860                         case *basicMatrix:
1861                                 err = (*Dense)(aU).Inverse(test.a)
1862                         default:
1863                                 continue
1864                         }
1865                         m.Mul(aU, &tmp)
1866                 } else {
1867                         switch a := test.a.(type) {
1868                         case *Dense:
1869                                 err = a.Inverse(test.a)
1870                                 m.Mul(a, &tmp)
1871                         case *basicMatrix:
1872                                 err = (*Dense)(a).Inverse(test.a)
1873                                 m.Mul(a, &tmp)
1874                         default:
1875                                 continue
1876                         }
1877                 }
1878                 if err != nil {
1879                         t.Errorf("Error computing inverse: %v", err)
1880                 }
1881                 if !equalApprox(eye, &m, 1e-14, false) {
1882                         t.Errorf("Case %d, A^-1 * A != I", i)
1883                         fmt.Println(Formatted(&m))
1884                 }
1885         }
1886 }
1887
1888 var (
1889         wd *Dense
1890 )
1891
1892 func BenchmarkMulDense100Half(b *testing.B)        { denseMulBench(b, 100, 0.5) }
1893 func BenchmarkMulDense100Tenth(b *testing.B)       { denseMulBench(b, 100, 0.1) }
1894 func BenchmarkMulDense1000Half(b *testing.B)       { denseMulBench(b, 1000, 0.5) }
1895 func BenchmarkMulDense1000Tenth(b *testing.B)      { denseMulBench(b, 1000, 0.1) }
1896 func BenchmarkMulDense1000Hundredth(b *testing.B)  { denseMulBench(b, 1000, 0.01) }
1897 func BenchmarkMulDense1000Thousandth(b *testing.B) { denseMulBench(b, 1000, 0.001) }
1898 func denseMulBench(b *testing.B, size int, rho float64) {
1899         b.StopTimer()
1900         a, _ := randDense(size, rho, rand.NormFloat64)
1901         d, _ := randDense(size, rho, rand.NormFloat64)
1902         b.StartTimer()
1903         for i := 0; i < b.N; i++ {
1904                 var n Dense
1905                 n.Mul(a, d)
1906                 wd = &n
1907         }
1908 }
1909
1910 func BenchmarkPreMulDense100Half(b *testing.B)        { densePreMulBench(b, 100, 0.5) }
1911 func BenchmarkPreMulDense100Tenth(b *testing.B)       { densePreMulBench(b, 100, 0.1) }
1912 func BenchmarkPreMulDense1000Half(b *testing.B)       { densePreMulBench(b, 1000, 0.5) }
1913 func BenchmarkPreMulDense1000Tenth(b *testing.B)      { densePreMulBench(b, 1000, 0.1) }
1914 func BenchmarkPreMulDense1000Hundredth(b *testing.B)  { densePreMulBench(b, 1000, 0.01) }
1915 func BenchmarkPreMulDense1000Thousandth(b *testing.B) { densePreMulBench(b, 1000, 0.001) }
1916 func densePreMulBench(b *testing.B, size int, rho float64) {
1917         b.StopTimer()
1918         a, _ := randDense(size, rho, rand.NormFloat64)
1919         d, _ := randDense(size, rho, rand.NormFloat64)
1920         wd = NewDense(size, size, nil)
1921         b.StartTimer()
1922         for i := 0; i < b.N; i++ {
1923                 wd.Mul(a, d)
1924         }
1925 }
1926
1927 func BenchmarkRow10(b *testing.B)   { rowBench(b, 10) }
1928 func BenchmarkRow100(b *testing.B)  { rowBench(b, 100) }
1929 func BenchmarkRow1000(b *testing.B) { rowBench(b, 1000) }
1930
1931 func rowBench(b *testing.B, size int) {
1932         a, _ := randDense(size, 1, rand.NormFloat64)
1933         _, c := a.Dims()
1934         dst := make([]float64, c)
1935
1936         b.ResetTimer()
1937         for i := 0; i < b.N; i++ {
1938                 Row(dst, 0, a)
1939         }
1940 }
1941
1942 func BenchmarkExp10(b *testing.B)   { expBench(b, 10) }
1943 func BenchmarkExp100(b *testing.B)  { expBench(b, 100) }
1944 func BenchmarkExp1000(b *testing.B) { expBench(b, 1000) }
1945
1946 func expBench(b *testing.B, size int) {
1947         a, _ := randDense(size, 1, rand.NormFloat64)
1948
1949         b.ResetTimer()
1950         var m Dense
1951         for i := 0; i < b.N; i++ {
1952                 m.Exp(a)
1953         }
1954 }
1955
1956 func BenchmarkPow10_3(b *testing.B)   { powBench(b, 10, 3) }
1957 func BenchmarkPow100_3(b *testing.B)  { powBench(b, 100, 3) }
1958 func BenchmarkPow1000_3(b *testing.B) { powBench(b, 1000, 3) }
1959 func BenchmarkPow10_4(b *testing.B)   { powBench(b, 10, 4) }
1960 func BenchmarkPow100_4(b *testing.B)  { powBench(b, 100, 4) }
1961 func BenchmarkPow1000_4(b *testing.B) { powBench(b, 1000, 4) }
1962 func BenchmarkPow10_5(b *testing.B)   { powBench(b, 10, 5) }
1963 func BenchmarkPow100_5(b *testing.B)  { powBench(b, 100, 5) }
1964 func BenchmarkPow1000_5(b *testing.B) { powBench(b, 1000, 5) }
1965 func BenchmarkPow10_6(b *testing.B)   { powBench(b, 10, 6) }
1966 func BenchmarkPow100_6(b *testing.B)  { powBench(b, 100, 6) }
1967 func BenchmarkPow1000_6(b *testing.B) { powBench(b, 1000, 6) }
1968 func BenchmarkPow10_7(b *testing.B)   { powBench(b, 10, 7) }
1969 func BenchmarkPow100_7(b *testing.B)  { powBench(b, 100, 7) }
1970 func BenchmarkPow1000_7(b *testing.B) { powBench(b, 1000, 7) }
1971 func BenchmarkPow10_8(b *testing.B)   { powBench(b, 10, 8) }
1972 func BenchmarkPow100_8(b *testing.B)  { powBench(b, 100, 8) }
1973 func BenchmarkPow1000_8(b *testing.B) { powBench(b, 1000, 8) }
1974 func BenchmarkPow10_9(b *testing.B)   { powBench(b, 10, 9) }
1975 func BenchmarkPow100_9(b *testing.B)  { powBench(b, 100, 9) }
1976 func BenchmarkPow1000_9(b *testing.B) { powBench(b, 1000, 9) }
1977
1978 func powBench(b *testing.B, size, n int) {
1979         a, _ := randDense(size, 1, rand.NormFloat64)
1980
1981         b.ResetTimer()
1982         var m Dense
1983         for i := 0; i < b.N; i++ {
1984                 m.Pow(a, n)
1985         }
1986 }
1987
1988 func BenchmarkMulTransDense100Half(b *testing.B)        { denseMulTransBench(b, 100, 0.5) }
1989 func BenchmarkMulTransDense100Tenth(b *testing.B)       { denseMulTransBench(b, 100, 0.1) }
1990 func BenchmarkMulTransDense1000Half(b *testing.B)       { denseMulTransBench(b, 1000, 0.5) }
1991 func BenchmarkMulTransDense1000Tenth(b *testing.B)      { denseMulTransBench(b, 1000, 0.1) }
1992 func BenchmarkMulTransDense1000Hundredth(b *testing.B)  { denseMulTransBench(b, 1000, 0.01) }
1993 func BenchmarkMulTransDense1000Thousandth(b *testing.B) { denseMulTransBench(b, 1000, 0.001) }
1994 func denseMulTransBench(b *testing.B, size int, rho float64) {
1995         b.StopTimer()
1996         a, _ := randDense(size, rho, rand.NormFloat64)
1997         d, _ := randDense(size, rho, rand.NormFloat64)
1998         b.StartTimer()
1999         for i := 0; i < b.N; i++ {
2000                 var n Dense
2001                 n.Mul(a, d.T())
2002                 wd = &n
2003         }
2004 }
2005
2006 func BenchmarkMulTransDenseSym100Half(b *testing.B)        { denseMulTransSymBench(b, 100, 0.5) }
2007 func BenchmarkMulTransDenseSym100Tenth(b *testing.B)       { denseMulTransSymBench(b, 100, 0.1) }
2008 func BenchmarkMulTransDenseSym1000Half(b *testing.B)       { denseMulTransSymBench(b, 1000, 0.5) }
2009 func BenchmarkMulTransDenseSym1000Tenth(b *testing.B)      { denseMulTransSymBench(b, 1000, 0.1) }
2010 func BenchmarkMulTransDenseSym1000Hundredth(b *testing.B)  { denseMulTransSymBench(b, 1000, 0.01) }
2011 func BenchmarkMulTransDenseSym1000Thousandth(b *testing.B) { denseMulTransSymBench(b, 1000, 0.001) }
2012 func denseMulTransSymBench(b *testing.B, size int, rho float64) {
2013         b.StopTimer()
2014         a, _ := randDense(size, rho, rand.NormFloat64)
2015         b.StartTimer()
2016         for i := 0; i < b.N; i++ {
2017                 var n Dense
2018                 n.Mul(a, a.T())
2019                 wd = &n
2020         }
2021 }