github.com/gopherd/gonum@v0.0.4/mat/symmetric_test.go (about)

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