diff --git a/mat/cholesky_test.go b/mat/cholesky_test.go
index a8325045..0b753181 100644
--- a/mat/cholesky_test.go
+++ b/mat/cholesky_test.go
@@ -331,10 +331,11 @@ func TestCloneCholesky(t *testing.T) {
 
 func TestCholeskyInverseTo(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, n := range []int{1, 3, 5, 9} {
 		data := make([]float64, n*n)
 		for i := range data {
-			data[i] = rand.NormFloat64()
+			data[i] = rnd.NormFloat64()
 		}
 		var s SymDense
 		s.SymOuterK(1, NewDense(n, n, data))
@@ -363,13 +364,13 @@ func TestCholeskyInverseTo(t *testing.T) {
 
 func TestCholeskySymRankOne(t *testing.T) {
 	t.Parallel()
-	rand.Seed(1)
+	rnd := rand.New(rand.NewSource(1))
 	for _, n := range []int{1, 2, 3, 4, 5, 7, 10, 20, 50, 100} {
 		for k := 0; k < 50; k++ {
 			// Construct a random positive definite matrix.
 			data := make([]float64, n*n)
 			for i := range data {
-				data[i] = rand.NormFloat64()
+				data[i] = rnd.NormFloat64()
 			}
 			var a SymDense
 			a.SymOuterK(1, NewDense(n, n, data))
@@ -377,10 +378,10 @@ func TestCholeskySymRankOne(t *testing.T) {
 			// Construct random data for updating.
 			xdata := make([]float64, n)
 			for i := range xdata {
-				xdata[i] = rand.NormFloat64()
+				xdata[i] = rnd.NormFloat64()
 			}
 			x := NewVecDense(n, xdata)
-			alpha := rand.NormFloat64()
+			alpha := rnd.NormFloat64()
 
 			// Compute the updated matrix directly. If alpha > 0, there are no
 			// issues. If alpha < 0, it could be that the final matrix is not
diff --git a/mat/dense_test.go b/mat/dense_test.go
index 46a3c11f..f6a32e02 100644
--- a/mat/dense_test.go
+++ b/mat/dense_test.go
@@ -937,7 +937,7 @@ func TestDenseMul(t *testing.T) {
 	testTwoInput(t, "Mul", &Dense{}, method, denseComparison, legalTypesAll, legalSizeMul, 1e-14)
 }
 
-func randDense(size int, rho float64, rnd func() float64) (*Dense, error) {
+func randDense(size int, rho float64, src rand.Source) (*Dense, error) {
 	if size == 0 {
 		return nil, ErrZeroLength
 	}
@@ -948,10 +948,11 @@ func randDense(size int, rho float64, rnd func() float64) (*Dense, error) {
 		},
 		capRows: size, capCols: size,
 	}
+	rnd := rand.New(src)
 	for i := 0; i < size; i++ {
 		for j := 0; j < size; j++ {
-			if rand.Float64() < rho {
-				d.Set(i, j, rnd())
+			if rnd.Float64() < rho {
+				d.Set(i, j, rnd.NormFloat64())
 			}
 		}
 	}
@@ -2099,9 +2100,10 @@ func BenchmarkMulDense1000Tenth(b *testing.B)      { denseMulBench(b, 1000, 0.1)
 func BenchmarkMulDense1000Hundredth(b *testing.B)  { denseMulBench(b, 1000, 0.01) }
 func BenchmarkMulDense1000Thousandth(b *testing.B) { denseMulBench(b, 1000, 0.001) }
 func denseMulBench(b *testing.B, size int, rho float64) {
+	src := rand.NewSource(1)
 	b.StopTimer()
-	a, _ := randDense(size, rho, rand.NormFloat64)
-	d, _ := randDense(size, rho, rand.NormFloat64)
+	a, _ := randDense(size, rho, src)
+	d, _ := randDense(size, rho, src)
 	b.StartTimer()
 	for i := 0; i < b.N; i++ {
 		var n Dense
@@ -2117,9 +2119,10 @@ func BenchmarkPreMulDense1000Tenth(b *testing.B)      { densePreMulBench(b, 1000
 func BenchmarkPreMulDense1000Hundredth(b *testing.B)  { densePreMulBench(b, 1000, 0.01) }
 func BenchmarkPreMulDense1000Thousandth(b *testing.B) { densePreMulBench(b, 1000, 0.001) }
 func densePreMulBench(b *testing.B, size int, rho float64) {
+	src := rand.NewSource(1)
 	b.StopTimer()
-	a, _ := randDense(size, rho, rand.NormFloat64)
-	d, _ := randDense(size, rho, rand.NormFloat64)
+	a, _ := randDense(size, rho, src)
+	d, _ := randDense(size, rho, src)
 	wd = NewDense(size, size, nil)
 	b.StartTimer()
 	for i := 0; i < b.N; i++ {
@@ -2132,7 +2135,8 @@ func BenchmarkDenseRow100(b *testing.B)  { rowDenseBench(b, 100) }
 func BenchmarkDenseRow1000(b *testing.B) { rowDenseBench(b, 1000) }
 
 func rowDenseBench(b *testing.B, size int) {
-	a, _ := randDense(size, 1, rand.NormFloat64)
+	src := rand.NewSource(1)
+	a, _ := randDense(size, 1, src)
 	_, c := a.Dims()
 	dst := make([]float64, c)
 
@@ -2147,7 +2151,8 @@ func BenchmarkDenseExp100(b *testing.B)  { expDenseBench(b, 100) }
 func BenchmarkDenseExp1000(b *testing.B) { expDenseBench(b, 1000) }
 
 func expDenseBench(b *testing.B, size int) {
-	a, _ := randDense(size, 1, rand.NormFloat64)
+	src := rand.NewSource(1)
+	a, _ := randDense(size, 1, src)
 
 	b.ResetTimer()
 	var m Dense
@@ -2179,7 +2184,8 @@ func BenchmarkDensePow100_9(b *testing.B)  { powDenseBench(b, 100, 9) }
 func BenchmarkDensePow1000_9(b *testing.B) { powDenseBench(b, 1000, 9) }
 
 func powDenseBench(b *testing.B, size, n int) {
-	a, _ := randDense(size, 1, rand.NormFloat64)
+	src := rand.NewSource(1)
+	a, _ := randDense(size, 1, src)
 
 	b.ResetTimer()
 	var m Dense
@@ -2195,9 +2201,10 @@ func BenchmarkDenseMulTransDense1000Tenth(b *testing.B)      { denseMulTransBenc
 func BenchmarkDenseMulTransDense1000Hundredth(b *testing.B)  { denseMulTransBench(b, 1000, 0.01) }
 func BenchmarkDenseMulTransDense1000Thousandth(b *testing.B) { denseMulTransBench(b, 1000, 0.001) }
 func denseMulTransBench(b *testing.B, size int, rho float64) {
+	src := rand.NewSource(1)
 	b.StopTimer()
-	a, _ := randDense(size, rho, rand.NormFloat64)
-	d, _ := randDense(size, rho, rand.NormFloat64)
+	a, _ := randDense(size, rho, src)
+	d, _ := randDense(size, rho, src)
 	b.StartTimer()
 	for i := 0; i < b.N; i++ {
 		var n Dense
@@ -2215,8 +2222,9 @@ func BenchmarkDenseMulTransDenseSym1000Thousandth(b *testing.B) {
 	denseMulTransSymBench(b, 1000, 0.001)
 }
 func denseMulTransSymBench(b *testing.B, size int, rho float64) {
+	src := rand.NewSource(1)
 	b.StopTimer()
-	a, _ := randDense(size, rho, rand.NormFloat64)
+	a, _ := randDense(size, rho, src)
 	b.StartTimer()
 	for i := 0; i < b.N; i++ {
 		var n Dense
@@ -2230,7 +2238,8 @@ func BenchmarkDenseSum1000(b *testing.B) { denseSumBench(b, 1000) }
 var denseSumForBench float64
 
 func denseSumBench(b *testing.B, size int) {
-	a, _ := randDense(size, 1.0, rand.NormFloat64)
+	src := rand.NewSource(1)
+	a, _ := randDense(size, 1.0, src)
 	b.ResetTimer()
 	for i := 0; i < b.N; i++ {
 		denseSumForBench = Sum(a)
diff --git a/mat/inner_test.go b/mat/inner_test.go
index ad53076e..269b2ea3 100644
--- a/mat/inner_test.go
+++ b/mat/inner_test.go
@@ -8,6 +8,8 @@ import (
 	"math"
 	"testing"
 
+	"golang.org/x/exp/rand"
+
 	"gonum.org/v1/gonum/blas/blas64"
 	"gonum.org/v1/gonum/blas/testblas"
 )
@@ -153,12 +155,13 @@ func makeVecDenseInc(inc int, f []float64) *VecDense {
 }
 
 func benchmarkInner(b *testing.B, m, n int) {
+	src := rand.NewSource(1)
 	x := NewVecDense(m, nil)
-	randomSlice(x.mat.Data)
+	randomSlice(x.mat.Data, src)
 	y := NewVecDense(n, nil)
-	randomSlice(y.mat.Data)
+	randomSlice(y.mat.Data, src)
 	data := make([]float64, m*n)
-	randomSlice(data)
+	randomSlice(data, src)
 	mat := &Dense{mat: blas64.General{Rows: m, Cols: n, Stride: n, Data: data}, capRows: m, capCols: n}
 	b.ResetTimer()
 	for i := 0; i < b.N; i++ {
diff --git a/mat/list_test.go b/mat/list_test.go
index 56e8f30f..24bb6191 100644
--- a/mat/list_test.go
+++ b/mat/list_test.go
@@ -342,13 +342,14 @@ func retranspose(a, m Matrix) Matrix {
 }
 
 // makeRandOf returns a new randomly filled m×n matrix of the underlying matrix type.
-func makeRandOf(a Matrix, m, n int) Matrix {
+func makeRandOf(a Matrix, m, n int, src rand.Source) Matrix {
+	rnd := rand.New(src)
 	var rMatrix Matrix
 	switch t := a.(type) {
 	default:
 		panic("unknown type for make rand of")
 	case Untransposer:
-		rMatrix = retranspose(a, makeRandOf(t.Untranspose(), n, m))
+		rMatrix = retranspose(a, makeRandOf(t.Untranspose(), n, m, src))
 	case *Dense, *basicMatrix:
 		var mat = &Dense{}
 		if m != 0 && n != 0 {
@@ -356,7 +357,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 		}
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
-				mat.Set(i, j, rand.NormFloat64())
+				mat.Set(i, j, rnd.NormFloat64())
 			}
 		}
 		rMatrix = returnAs(mat, t)
@@ -380,7 +381,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 			},
 		}
 		for i := 0; i < length; i++ {
-			mat.SetVec(i, rand.NormFloat64())
+			mat.SetVec(i, rnd.NormFloat64())
 		}
 		return mat
 	case *basicVector:
@@ -394,7 +395,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 			m: make([]float64, m),
 		}
 		for i := 0; i < m; i++ {
-			mat.m[i] = rand.NormFloat64()
+			mat.m[i] = rnd.NormFloat64()
 		}
 		return mat
 	case *SymDense, *basicSymmetric:
@@ -407,7 +408,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 		}
 		for i := 0; i < m; i++ {
 			for j := i; j < n; j++ {
-				mat.SetSym(i, j, rand.NormFloat64())
+				mat.SetSym(i, j, rnd.NormFloat64())
 			}
 		}
 		rMatrix = returnAs(mat, t)
@@ -439,13 +440,13 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 		if triKind == Upper {
 			for i := 0; i < m; i++ {
 				for j := i; j < n; j++ {
-					mat.SetTri(i, j, rand.NormFloat64())
+					mat.SetTri(i, j, rnd.NormFloat64())
 				}
 			}
 		} else {
 			for i := 0; i < m; i++ {
 				for j := 0; j <= i; j++ {
-					mat.SetTri(i, j, rand.NormFloat64())
+					mat.SetTri(i, j, rnd.NormFloat64())
 				}
 			}
 		}
@@ -464,7 +465,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 		kl = min(kl, m-1)
 		data := make([]float64, min(m, n+kl)*(kl+ku+1))
 		for i := range data {
-			data[i] = rand.NormFloat64()
+			data[i] = rnd.NormFloat64()
 		}
 		mat := NewBandDense(m, n, kl, ku, data)
 		rMatrix = returnAs(mat, t)
@@ -482,7 +483,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 		k = min(k, m-1) // Special case for small sizes.
 		data := make([]float64, m*(k+1))
 		for i := range data {
-			data[i] = rand.NormFloat64()
+			data[i] = rnd.NormFloat64()
 		}
 		mat := NewSymBandDense(n, k, data)
 		rMatrix = returnAs(mat, t)
@@ -503,7 +504,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 		k = min(k, m-1) // Special case for small sizes.
 		data := make([]float64, m*(k+1))
 		for i := range data {
-			data[i] = rand.NormFloat64()
+			data[i] = rnd.NormFloat64()
 		}
 		mat := NewTriBandDense(n, k, triKind, data)
 		rMatrix = returnAs(mat, t)
@@ -529,7 +530,7 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 			},
 		}
 		for i := 0; i < n; i++ {
-			mat.SetDiag(i, rand.Float64())
+			mat.SetDiag(i, rnd.Float64())
 		}
 		rMatrix = returnAs(mat, t)
 	}
@@ -908,6 +909,7 @@ func testOneInputFunc(t *testing.T,
 	// legalSize returns true if the size is valid for the function.
 	legalSize func(r, c int) bool,
 ) {
+	src := rand.NewSource(1)
 	for _, aMat := range testMatrices {
 		for _, test := range sizes {
 			// Skip the test if the argument would not be assignable to the
@@ -919,7 +921,7 @@ func testOneInputFunc(t *testing.T,
 			if !legalDims(aMat, test.ar, test.ac) {
 				continue
 			}
-			a := makeRandOf(aMat, test.ar, test.ac)
+			a := makeRandOf(aMat, test.ar, test.ac, src)
 
 			// Compute the true answer if the sizes are legal.
 			dimsOK := legalSize(test.ar, test.ac)
@@ -1072,6 +1074,7 @@ func testTwoInputFunc(t *testing.T,
 	// legalSize returns true if the sizes are valid for the function.
 	legalSize func(ar, ac, br, bc int) bool,
 ) {
+	src := rand.NewSource(1)
 	for _, aMat := range testMatrices {
 		for _, bMat := range testMatrices {
 			// Loop over all of the size combinations (bigger, smaller, etc.).
@@ -1088,8 +1091,8 @@ func testTwoInputFunc(t *testing.T,
 				if !legalDims(bMat, test.br, test.bc) {
 					continue
 				}
-				a := makeRandOf(aMat, test.ar, test.ac)
-				b := makeRandOf(bMat, test.br, test.bc)
+				a := makeRandOf(aMat, test.ar, test.ac, src)
+				b := makeRandOf(bMat, test.br, test.bc, src)
 
 				// Compute the true answer if the sizes are legal.
 				dimsOK := legalSize(test.ar, test.ac, test.br, test.bc)
@@ -1158,6 +1161,7 @@ func testOneInput(t *testing.T,
 	// tol is the tolerance for equality when comparing method results.
 	tol float64,
 ) {
+	src := rand.NewSource(1)
 	for _, aMat := range testMatrices {
 		for _, test := range sizes {
 			// Skip the test if the argument would not be assignable to the
@@ -1169,7 +1173,7 @@ func testOneInput(t *testing.T,
 			if !legalDims(aMat, test.ar, test.ac) {
 				continue
 			}
-			a := makeRandOf(aMat, test.ar, test.ac)
+			a := makeRandOf(aMat, test.ar, test.ac, src)
 
 			// Compute the true answer if the sizes are legal.
 			dimsOK := legalSize(test.ar, test.ac)
@@ -1184,7 +1188,7 @@ func testOneInput(t *testing.T,
 			// Test the method for a zero-value of the receiver.
 			aType, aTrans := untranspose(a)
 			errStr := fmt.Sprintf("%T.%s(%T), size: %#v, atrans %v", receiver, name, aType, test, aTrans)
-			zero := makeRandOf(receiver, 0, 0)
+			zero := makeRandOf(receiver, 0, 0, src)
 			panicked, err := panics(func() { method(zero, a) })
 			if !dimsOK && !panicked {
 				t.Errorf("Did not panic with illegal size: %s", errStr)
@@ -1209,7 +1213,7 @@ func testOneInput(t *testing.T,
 			// The receiver has been overwritten in place so use its size
 			// to construct a new random matrix.
 			rr, rc := zero.Dims()
-			neverZero := makeRandOf(receiver, rr, rc)
+			neverZero := makeRandOf(receiver, rr, rc, src)
 			panicked, _ = panics(func() { method(neverZero, a) })
 			if panicked {
 				t.Errorf("Panicked with non-zero receiver: %s", errStr)
@@ -1223,26 +1227,26 @@ func testOneInput(t *testing.T,
 			default:
 				panic("matrix type not coded for incorrect receiver size")
 			case *Dense:
-				wrongSize := makeRandOf(receiver, rr+1, rc)
+				wrongSize := makeRandOf(receiver, rr+1, rc, src)
 				panicked, _ = panics(func() { method(wrongSize, a) })
 				if !panicked {
 					t.Errorf("Did not panic with wrong number of rows: %s", errStr)
 				}
-				wrongSize = makeRandOf(receiver, rr, rc+1)
+				wrongSize = makeRandOf(receiver, rr, rc+1, src)
 				panicked, _ = panics(func() { method(wrongSize, a) })
 				if !panicked {
 					t.Errorf("Did not panic with wrong number of columns: %s", errStr)
 				}
 			case *TriDense, *SymDense:
 				// Add to the square size.
-				wrongSize := makeRandOf(receiver, rr+1, rc+1)
+				wrongSize := makeRandOf(receiver, rr+1, rc+1, src)
 				panicked, _ = panics(func() { method(wrongSize, a) })
 				if !panicked {
 					t.Errorf("Did not panic with wrong size: %s", errStr)
 				}
 			case *VecDense:
 				// Add to the column length.
-				wrongSize := makeRandOf(receiver, rr+1, rc)
+				wrongSize := makeRandOf(receiver, rr+1, rc, src)
 				panicked, _ = panics(func() { method(wrongSize, a) })
 				if !panicked {
 					t.Errorf("Did not panic with wrong number of rows: %s", errStr)
@@ -1304,6 +1308,7 @@ func testTwoInput(t *testing.T,
 	// tol is the tolerance for equality when comparing method results.
 	tol float64,
 ) {
+	src := rand.NewSource(1)
 	for _, aMat := range testMatrices {
 		for _, bMat := range testMatrices {
 			// Loop over all of the size combinations (bigger, smaller, etc.).
@@ -1320,8 +1325,8 @@ func testTwoInput(t *testing.T,
 				if !legalDims(bMat, test.br, test.bc) {
 					continue
 				}
-				a := makeRandOf(aMat, test.ar, test.ac)
-				b := makeRandOf(bMat, test.br, test.bc)
+				a := makeRandOf(aMat, test.ar, test.ac, src)
+				b := makeRandOf(bMat, test.br, test.bc, src)
 
 				// Compute the true answer if the sizes are legal.
 				dimsOK := legalSize(test.ar, test.ac, test.br, test.bc)
@@ -1339,7 +1344,7 @@ func testTwoInput(t *testing.T,
 				aType, aTrans := untranspose(a)
 				bType, bTrans := untranspose(b)
 				errStr := fmt.Sprintf("%T.%s(%T, %T), sizes: %#v, atrans %v, btrans %v", receiver, name, aType, bType, test, aTrans, bTrans)
-				zero := makeRandOf(receiver, 0, 0)
+				zero := makeRandOf(receiver, 0, 0, src)
 				panicked, err := panics(func() { method(zero, a, b) })
 				if !dimsOK && !panicked {
 					t.Errorf("Did not panic with illegal size: %s", errStr)
@@ -1369,7 +1374,7 @@ func testTwoInput(t *testing.T,
 				// The receiver has been overwritten in place so use its size
 				// to construct a new random matrix.
 				rr, rc := wasZero.Dims()
-				neverZero := makeRandOf(receiver, rr, rc)
+				neverZero := makeRandOf(receiver, rr, rc, src)
 				panicked, message := panics(func() { method(neverZero, a, b) })
 				if panicked {
 					t.Errorf("Panicked with non-zero receiver: %s: %s", errStr, message)
@@ -1384,26 +1389,26 @@ func testTwoInput(t *testing.T,
 				default:
 					panic("matrix type not coded for incorrect receiver size")
 				case *Dense:
-					wrongSize := makeRandOf(receiver, rr+1, rc)
+					wrongSize := makeRandOf(receiver, rr+1, rc, src)
 					panicked, _ = panics(func() { method(wrongSize, a, b) })
 					if !panicked {
 						t.Errorf("Did not panic with wrong number of rows: %s", errStr)
 					}
-					wrongSize = makeRandOf(receiver, rr, rc+1)
+					wrongSize = makeRandOf(receiver, rr, rc+1, src)
 					panicked, _ = panics(func() { method(wrongSize, a, b) })
 					if !panicked {
 						t.Errorf("Did not panic with wrong number of columns: %s", errStr)
 					}
 				case *TriDense, *SymDense:
 					// Add to the square size.
-					wrongSize := makeRandOf(receiver, rr+1, rc+1)
+					wrongSize := makeRandOf(receiver, rr+1, rc+1, src)
 					panicked, _ = panics(func() { method(wrongSize, a, b) })
 					if !panicked {
 						t.Errorf("Did not panic with wrong size: %s", errStr)
 					}
 				case *VecDense:
 					// Add to the column length.
-					wrongSize := makeRandOf(receiver, rr+1, rc)
+					wrongSize := makeRandOf(receiver, rr+1, rc, src)
 					panicked, _ = panics(func() { method(wrongSize, a, b) })
 					if !panicked {
 						t.Errorf("Did not panic with wrong number of rows: %s", errStr)
@@ -1475,7 +1480,7 @@ func testTwoInput(t *testing.T,
 					// Compute the real answer for this case. It is different
 					// from the initial answer since now a and b have the
 					// same data.
-					zero = makeRandOf(wasZero, 0, 0)
+					zero = makeRandOf(wasZero, 0, 0, src)
 					method(zero, aSame, bSame)
 					wasZero, zero = zero, nil // Nil-out zero so we detect illegal use.
 					preData := underlyingData(receiver)
diff --git a/mat/lq_test.go b/mat/lq_test.go
index 3fcc6ef0..64dab97d 100644
--- a/mat/lq_test.go
+++ b/mat/lq_test.go
@@ -12,6 +12,7 @@ import (
 
 func TestLQ(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		m, n int
 	}{
@@ -23,7 +24,7 @@ func TestLQ(t *testing.T) {
 		a := NewDense(m, n, nil)
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
-				a.Set(i, j, rand.NormFloat64())
+				a.Set(i, j, rnd.NormFloat64())
 			}
 		}
 		var want Dense
@@ -50,6 +51,7 @@ func TestLQ(t *testing.T) {
 
 func TestLQSolveTo(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n, bc int
@@ -65,7 +67,7 @@ func TestLQSolveTo(t *testing.T) {
 			a := NewDense(m, n, nil)
 			for i := 0; i < m; i++ {
 				for j := 0; j < n; j++ {
-					a.Set(i, j, rand.Float64())
+					a.Set(i, j, rnd.Float64())
 				}
 			}
 			br := m
@@ -75,7 +77,7 @@ func TestLQSolveTo(t *testing.T) {
 			b := NewDense(br, bc, nil)
 			for i := 0; i < br; i++ {
 				for j := 0; j < bc; j++ {
-					b.Set(i, j, rand.Float64())
+					b.Set(i, j, rnd.Float64())
 				}
 			}
 			var x Dense
@@ -112,6 +114,7 @@ func TestLQSolveTo(t *testing.T) {
 
 func TestLQSolveToVec(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n int
@@ -124,7 +127,7 @@ func TestLQSolveToVec(t *testing.T) {
 			a := NewDense(m, n, nil)
 			for i := 0; i < m; i++ {
 				for j := 0; j < n; j++ {
-					a.Set(i, j, rand.Float64())
+					a.Set(i, j, rnd.Float64())
 				}
 			}
 			br := m
@@ -133,7 +136,7 @@ func TestLQSolveToVec(t *testing.T) {
 			}
 			b := NewVecDense(br, nil)
 			for i := 0; i < br; i++ {
-				b.SetVec(i, rand.Float64())
+				b.SetVec(i, rnd.Float64())
 			}
 			var x VecDense
 			lq := &LQ{}
diff --git a/mat/lu_test.go b/mat/lu_test.go
index 602b4056..e3fea10a 100644
--- a/mat/lu_test.go
+++ b/mat/lu_test.go
@@ -12,11 +12,12 @@ import (
 
 func TestLUD(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, n := range []int{1, 5, 10, 11, 50} {
 		a := NewDense(n, n, nil)
 		for i := 0; i < n; i++ {
 			for j := 0; j < n; j++ {
-				a.Set(i, j, rand.NormFloat64())
+				a.Set(i, j, rnd.NormFloat64())
 			}
 		}
 		var want Dense
@@ -41,6 +42,7 @@ func TestLUD(t *testing.T) {
 
 func TestLURankOne(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, pivoting := range []bool{true} {
 		for _, n := range []int{3, 10, 50} {
 			// Construct a random LU factorization
@@ -48,7 +50,7 @@ func TestLURankOne(t *testing.T) {
 			lu.lu = NewDense(n, n, nil)
 			for i := 0; i < n; i++ {
 				for j := 0; j < n; j++ {
-					lu.lu.Set(i, j, rand.Float64())
+					lu.lu.Set(i, j, rnd.Float64())
 				}
 			}
 			lu.pivot = make([]int, n)
@@ -59,18 +61,18 @@ func TestLURankOne(t *testing.T) {
 				// For each row, randomly swap with itself or a row after (like is done)
 				// in the actual LU factorization.
 				for i := range lu.pivot {
-					idx := i + rand.Intn(n-i)
+					idx := i + rnd.Intn(n-i)
 					lu.pivot[i], lu.pivot[idx] = lu.pivot[idx], lu.pivot[i]
 				}
 			}
 			// Apply a rank one update. Ensure the update magnitude is larger than
 			// the equal tolerance.
-			alpha := rand.Float64() + 1
+			alpha := rnd.Float64() + 1
 			x := NewVecDense(n, nil)
 			y := NewVecDense(n, nil)
 			for i := 0; i < n; i++ {
-				x.setVec(i, rand.Float64()+1)
-				y.setVec(i, rand.Float64()+1)
+				x.setVec(i, rnd.Float64()+1)
+				y.setVec(i, rnd.Float64()+1)
 			}
 			a := luReconstruct(lu)
 			a.RankOne(a, alpha, x, y)
@@ -109,6 +111,7 @@ func luReconstruct(lu *LU) *Dense {
 
 func TestLUSolveTo(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		n, bc int
 	}{
@@ -121,13 +124,13 @@ func TestLUSolveTo(t *testing.T) {
 		a := NewDense(n, n, nil)
 		for i := 0; i < n; i++ {
 			for j := 0; j < n; j++ {
-				a.Set(i, j, rand.NormFloat64())
+				a.Set(i, j, rnd.NormFloat64())
 			}
 		}
 		b := NewDense(n, bc, nil)
 		for i := 0; i < n; i++ {
 			for j := 0; j < bc; j++ {
-				b.Set(i, j, rand.NormFloat64())
+				b.Set(i, j, rnd.NormFloat64())
 			}
 		}
 		var lu LU
@@ -169,16 +172,17 @@ func TestLUSolveToCond(t *testing.T) {
 
 func TestLUSolveVecTo(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, n := range []int{5, 10} {
 		a := NewDense(n, n, nil)
 		for i := 0; i < n; i++ {
 			for j := 0; j < n; j++ {
-				a.Set(i, j, rand.NormFloat64())
+				a.Set(i, j, rnd.NormFloat64())
 			}
 		}
 		b := NewVecDense(n, nil)
 		for i := 0; i < n; i++ {
-			b.SetVec(i, rand.NormFloat64())
+			b.SetVec(i, rnd.NormFloat64())
 		}
 		var lu LU
 		lu.Factorize(a)
diff --git a/mat/mul_test.go b/mat/mul_test.go
index 5dafeac1..b15a49f7 100644
--- a/mat/mul_test.go
+++ b/mat/mul_test.go
@@ -17,6 +17,7 @@ import (
 // TODO: Need to add tests where one is overwritten.
 func TestMulTypes(t *testing.T) {
 	t.Parallel()
+	src := rand.NewSource(1)
 	for _, test := range []struct {
 		ar     int
 		ac     int
@@ -88,11 +89,11 @@ func TestMulTypes(t *testing.T) {
 
 		// Generate random matrices
 		avec := make([]float64, ar*ac)
-		randomSlice(avec)
+		randomSlice(avec, src)
 		a := NewDense(ar, ac, avec)
 
 		bvec := make([]float64, br*bc)
-		randomSlice(bvec)
+		randomSlice(bvec, src)
 
 		b := NewDense(br, bc, bvec)
 
@@ -132,7 +133,7 @@ func TestMulTypes(t *testing.T) {
 
 		// Normal multiply with existing receiver
 		c := NewDense(ar, bc, cvec)
-		randomSlice(cvec)
+		randomSlice(cvec, src)
 		testMul(t, a, b, c, acomp, bcomp, ccomp, false, "existing receiver")
 
 		// Cast a as a basic matrix
@@ -144,18 +145,19 @@ func TestMulTypes(t *testing.T) {
 		testMul(t, a, bm, d, acomp, bcomp, ccomp, true, "b is basic, receiver is empty")
 		d.Reset()
 		testMul(t, am, bm, d, acomp, bcomp, ccomp, true, "both basic, receiver is empty")
-		randomSlice(cvec)
+		randomSlice(cvec, src)
 		testMul(t, am, b, d, acomp, bcomp, ccomp, true, "a is basic, receiver is full")
-		randomSlice(cvec)
+		randomSlice(cvec, src)
 		testMul(t, a, bm, d, acomp, bcomp, ccomp, true, "b is basic, receiver is full")
-		randomSlice(cvec)
+		randomSlice(cvec, src)
 		testMul(t, am, bm, d, acomp, bcomp, ccomp, true, "both basic, receiver is full")
 	}
 }
 
-func randomSlice(s []float64) {
+func randomSlice(s []float64, src rand.Source) {
+	rnd := rand.New(src)
 	for i := range s {
-		s[i] = rand.NormFloat64()
+		s[i] = rnd.NormFloat64()
 	}
 }
 
diff --git a/mat/pool_test.go b/mat/pool_test.go
index 196d24ef..177cc1cc 100644
--- a/mat/pool_test.go
+++ b/mat/pool_test.go
@@ -73,9 +73,10 @@ func BenchmarkMulWorkspaceDense1000Tenth(b *testing.B)      { denseMulWorkspaceB
 func BenchmarkMulWorkspaceDense1000Hundredth(b *testing.B)  { denseMulWorkspaceBench(b, 1000, 0.01) }
 func BenchmarkMulWorkspaceDense1000Thousandth(b *testing.B) { denseMulWorkspaceBench(b, 1000, 0.001) }
 func denseMulWorkspaceBench(b *testing.B, size int, rho float64) {
+	src := rand.NewSource(1)
 	b.StopTimer()
-	a, _ := randDense(size, rho, rand.NormFloat64)
-	d, _ := randDense(size, rho, rand.NormFloat64)
+	a, _ := randDense(size, rho, src)
+	d, _ := randDense(size, rho, src)
 	b.StartTimer()
 	for i := 0; i < b.N; i++ {
 		a.Mul(a, d)
diff --git a/mat/product_test.go b/mat/product_test.go
index bf659099..97942ed9 100644
--- a/mat/product_test.go
+++ b/mat/product_test.go
@@ -92,6 +92,7 @@ var productTests = []struct {
 
 func TestProduct(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range productTests {
 		dimensions := test.factors
 		if dimensions == nil && test.n > 0 {
@@ -100,7 +101,7 @@ func TestProduct(t *testing.T) {
 				if i != 0 {
 					dimensions[i].r = dimensions[i-1].c
 				}
-				dimensions[i].c = rand.Intn(50) + 1
+				dimensions[i].c = rnd.Intn(50) + 1
 			}
 			dimensions[0].r = test.product.r
 			dimensions[test.n-1].c = test.product.c
@@ -109,7 +110,7 @@ func TestProduct(t *testing.T) {
 		for i, d := range dimensions {
 			data := make([]float64, d.r*d.c)
 			for i := range data {
-				data[i] = rand.Float64()
+				data[i] = rnd.Float64()
 			}
 			factors[i] = NewDense(d.r, d.c, data)
 		}
diff --git a/mat/qr_test.go b/mat/qr_test.go
index 4cbd173d..f3dfdc83 100644
--- a/mat/qr_test.go
+++ b/mat/qr_test.go
@@ -15,6 +15,7 @@ import (
 
 func TestQR(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		m, n int
 	}{
@@ -26,7 +27,7 @@ func TestQR(t *testing.T) {
 		a := NewDense(m, n, nil)
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
-				a.Set(i, j, rand.NormFloat64())
+				a.Set(i, j, rnd.NormFloat64())
 			}
 		}
 		var want Dense
@@ -74,6 +75,7 @@ func isOrthonormal(q *Dense, tol float64) bool {
 
 func TestQRSolveTo(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n, bc int
@@ -89,7 +91,7 @@ func TestQRSolveTo(t *testing.T) {
 			a := NewDense(m, n, nil)
 			for i := 0; i < m; i++ {
 				for j := 0; j < n; j++ {
-					a.Set(i, j, rand.Float64())
+					a.Set(i, j, rnd.Float64())
 				}
 			}
 			br := m
@@ -99,7 +101,7 @@ func TestQRSolveTo(t *testing.T) {
 			b := NewDense(br, bc, nil)
 			for i := 0; i < br; i++ {
 				for j := 0; j < bc; j++ {
-					b.Set(i, j, rand.Float64())
+					b.Set(i, j, rnd.Float64())
 				}
 			}
 			var x Dense
@@ -136,6 +138,7 @@ func TestQRSolveTo(t *testing.T) {
 
 func TestQRSolveVecTo(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n int
@@ -148,7 +151,7 @@ func TestQRSolveVecTo(t *testing.T) {
 			a := NewDense(m, n, nil)
 			for i := 0; i < m; i++ {
 				for j := 0; j < n; j++ {
-					a.Set(i, j, rand.Float64())
+					a.Set(i, j, rnd.Float64())
 				}
 			}
 			br := m
@@ -157,7 +160,7 @@ func TestQRSolveVecTo(t *testing.T) {
 			}
 			b := NewVecDense(br, nil)
 			for i := 0; i < br; i++ {
-				b.SetVec(i, rand.Float64())
+				b.SetVec(i, rnd.Float64())
 			}
 			var x VecDense
 			var qr QR
diff --git a/mat/solve_test.go b/mat/solve_test.go
index bb28532a..043b3021 100644
--- a/mat/solve_test.go
+++ b/mat/solve_test.go
@@ -12,6 +12,7 @@ import (
 
 func TestSolve(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	// Hand-coded cases.
 	for _, test := range []struct {
 		a         [][]float64
@@ -204,14 +205,14 @@ func TestSolve(t *testing.T) {
 		a := NewDense(m, n, nil)
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
-				a.Set(i, j, rand.Float64())
+				a.Set(i, j, rnd.Float64())
 			}
 		}
 		br := m
 		b := NewDense(br, bc, nil)
 		for i := 0; i < br; i++ {
 			for j := 0; j < bc; j++ {
-				b.Set(i, j, rand.Float64())
+				b.Set(i, j, rnd.Float64())
 			}
 		}
 		var x Dense
@@ -247,6 +248,7 @@ func TestSolve(t *testing.T) {
 
 func TestSolveVec(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		m, n int
 	}{
@@ -265,13 +267,13 @@ func TestSolveVec(t *testing.T) {
 		a := NewDense(m, n, nil)
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
-				a.Set(i, j, rand.Float64())
+				a.Set(i, j, rnd.Float64())
 			}
 		}
 		br := m
 		b := NewVecDense(br, nil)
 		for i := 0; i < br; i++ {
-			b.SetVec(i, rand.Float64())
+			b.SetVec(i, rnd.Float64())
 		}
 		var x VecDense
 		err := x.SolveVec(a, b)
diff --git a/mat/svd_test.go b/mat/svd_test.go
index c5d62dba..d40b6a65 100644
--- a/mat/svd_test.go
+++ b/mat/svd_test.go
@@ -14,6 +14,7 @@ import (
 
 func TestSVD(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	// Hand coded tests
 	for _, test := range []struct {
 		a *Dense
@@ -106,7 +107,7 @@ func TestSVD(t *testing.T) {
 		for trial := 0; trial < 10; trial++ {
 			a := NewDense(m, n, nil)
 			for i := range a.mat.Data {
-				a.mat.Data[i] = rand.NormFloat64()
+				a.mat.Data[i] = rnd.NormFloat64()
 			}
 			aCopy := DenseCopyOf(a)
 
diff --git a/mat/symmetric_test.go b/mat/symmetric_test.go
index 867e47e3..b8252b32 100644
--- a/mat/symmetric_test.go
+++ b/mat/symmetric_test.go
@@ -177,6 +177,7 @@ func TestSymDiagView(t *testing.T) {
 
 func TestSymAdd(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		n int
 	}{
@@ -190,11 +191,11 @@ func TestSymAdd(t *testing.T) {
 		n := test.n
 		a := NewSymDense(n, nil)
 		for i := range a.mat.Data {
-			a.mat.Data[i] = rand.Float64()
+			a.mat.Data[i] = rnd.Float64()
 		}
 		b := NewSymDense(n, nil)
 		for i := range a.mat.Data {
-			b.mat.Data[i] = rand.Float64()
+			b.mat.Data[i] = rnd.Float64()
 		}
 		var m Dense
 		m.Add(a, b)
@@ -239,6 +240,7 @@ func TestSymAdd(t *testing.T) {
 
 func TestCopy(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		n int
 	}{
@@ -252,7 +254,7 @@ func TestCopy(t *testing.T) {
 		n := test.n
 		a := NewSymDense(n, nil)
 		for i := range a.mat.Data {
-			a.mat.Data[i] = rand.Float64()
+			a.mat.Data[i] = rnd.Float64()
 		}
 		s := NewSymDense(n, nil)
 		s.CopySym(a)
@@ -287,6 +289,7 @@ func TestSymCopyPanic(t *testing.T) {
 
 func TestSymRankOne(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	const tol = 1e-15
 
 	for _, test := range []struct {
@@ -303,11 +306,11 @@ func TestSymRankOne(t *testing.T) {
 		alpha := 2.0
 		a := NewSymDense(n, nil)
 		for i := range a.mat.Data {
-			a.mat.Data[i] = rand.Float64()
+			a.mat.Data[i] = rnd.Float64()
 		}
 		x := make([]float64, n)
 		for i := range x {
-			x[i] = rand.Float64()
+			x[i] = rnd.Float64()
 		}
 
 		xMat := NewDense(n, 1, x)
@@ -386,6 +389,7 @@ func TestIssue250SymRankOne(t *testing.T) {
 
 func TestRankTwo(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		n int
 	}{
@@ -400,13 +404,13 @@ func TestRankTwo(t *testing.T) {
 		alpha := 2.0
 		a := NewSymDense(n, nil)
 		for i := range a.mat.Data {
-			a.mat.Data[i] = rand.Float64()
+			a.mat.Data[i] = rnd.Float64()
 		}
 		x := make([]float64, n)
 		y := make([]float64, n)
 		for i := range x {
-			x[i] = rand.Float64()
-			y[i] = rand.Float64()
+			x[i] = rnd.Float64()
+			y[i] = rnd.Float64()
 		}
 
 		xMat := NewDense(n, 1, x)
@@ -763,19 +767,21 @@ func BenchmarkSymSum1000(b *testing.B) { symSumBench(b, 1000) }
 var symSumForBench float64
 
 func symSumBench(b *testing.B, size int) {
-	a := randSymDense(size)
+	src := rand.NewSource(1)
+	a := randSymDense(size, src)
 	b.ResetTimer()
 	for i := 0; i < b.N; i++ {
 		symSumForBench = Sum(a)
 	}
 }
 
-func randSymDense(size int) *SymDense {
+func randSymDense(size int, src rand.Source) *SymDense {
+	rnd := rand.New(src)
 	backData := make([]float64, size*size)
 	for i := 0; i < size; i++ {
-		backData[i*size+i] = rand.Float64()
+		backData[i*size+i] = rnd.Float64()
 		for j := i + 1; j < size; j++ {
-			v := rand.Float64()
+			v := rnd.Float64()
 			backData[i*size+j] = v
 			backData[j*size+i] = v
 		}
diff --git a/mat/triangular_test.go b/mat/triangular_test.go
index 191babf0..2aca32a1 100644
--- a/mat/triangular_test.go
+++ b/mat/triangular_test.go
@@ -204,9 +204,11 @@ func TestTriDiagView(t *testing.T) {
 
 func TestTriDenseCopy(t *testing.T) {
 	t.Parallel()
+	src := rand.NewSource(1)
+	rnd := rand.New(src)
 	for i := 0; i < 100; i++ {
-		size := rand.Intn(100)
-		r, err := randDense(size, 0.9, rand.NormFloat64)
+		size := rnd.Intn(100)
+		r, err := randDense(size, 0.9, src)
 		if size == 0 {
 			if err != ErrZeroLength {
 				t.Fatalf("expected error %v: got: %v", ErrZeroLength, err)
@@ -255,9 +257,11 @@ func TestTriDenseCopy(t *testing.T) {
 
 func TestTriTriDenseCopy(t *testing.T) {
 	t.Parallel()
+	src := rand.NewSource(1)
+	rnd := rand.New(src)
 	for i := 0; i < 100; i++ {
-		size := rand.Intn(100)
-		r, err := randDense(size, 1, rand.NormFloat64)
+		size := rnd.Intn(100)
+		r, err := randDense(size, 1, src)
 		if size == 0 {
 			if err != ErrZeroLength {
 				t.Fatalf("expected error %v: got: %v", ErrZeroLength, err)
@@ -308,11 +312,12 @@ func TestTriTriDenseCopy(t *testing.T) {
 
 func TestTriInverse(t *testing.T) {
 	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for _, kind := range []TriKind{Upper, Lower} {
 		for _, n := range []int{1, 3, 5, 9} {
 			data := make([]float64, n*n)
 			for i := range data {
-				data[i] = rand.NormFloat64()
+				data[i] = rnd.NormFloat64()
 			}
 			a := NewTriDense(n, kind, data)
 			var tr TriDense
diff --git a/mat/vector_test.go b/mat/vector_test.go
index 75842806..4c54ac21 100644
--- a/mat/vector_test.go
+++ b/mat/vector_test.go
@@ -498,8 +498,9 @@ func BenchmarkAddScaledVec1000Inc20(b *testing.B)   { addScaledVecBench(b, 1000,
 func BenchmarkAddScaledVec10000Inc20(b *testing.B)  { addScaledVecBench(b, 10000, 20) }
 func BenchmarkAddScaledVec100000Inc20(b *testing.B) { addScaledVecBench(b, 100000, 20) }
 func addScaledVecBench(b *testing.B, size, inc int) {
-	x := randVecDense(size, inc, 1, rand.NormFloat64)
-	y := randVecDense(size, inc, 1, rand.NormFloat64)
+	src := rand.NewSource(1)
+	x := randVecDense(size, inc, 1, src)
+	y := randVecDense(size, inc, 1, src)
 	b.ResetTimer()
 	var v VecDense
 	for i := 0; i < b.N; i++ {
@@ -523,7 +524,8 @@ func BenchmarkScaleVec1000Inc20(b *testing.B)   { scaleVecBench(b, 1000, 20) }
 func BenchmarkScaleVec10000Inc20(b *testing.B)  { scaleVecBench(b, 10000, 20) }
 func BenchmarkScaleVec100000Inc20(b *testing.B) { scaleVecBench(b, 100000, 20) }
 func scaleVecBench(b *testing.B, size, inc int) {
-	x := randVecDense(size, inc, 1, rand.NormFloat64)
+	src := rand.NewSource(1)
+	x := randVecDense(size, inc, 1, src)
 	b.ResetTimer()
 	var v VecDense
 	for i := 0; i < b.N; i++ {
@@ -547,8 +549,9 @@ func BenchmarkAddVec1000Inc20(b *testing.B)   { addVecBench(b, 1000, 20) }
 func BenchmarkAddVec10000Inc20(b *testing.B)  { addVecBench(b, 10000, 20) }
 func BenchmarkAddVec100000Inc20(b *testing.B) { addVecBench(b, 100000, 20) }
 func addVecBench(b *testing.B, size, inc int) {
-	x := randVecDense(size, inc, 1, rand.NormFloat64)
-	y := randVecDense(size, inc, 1, rand.NormFloat64)
+	src := rand.NewSource(1)
+	x := randVecDense(size, inc, 1, src)
+	y := randVecDense(size, inc, 1, src)
 	b.ResetTimer()
 	var v VecDense
 	for i := 0; i < b.N; i++ {
@@ -572,8 +575,9 @@ func BenchmarkSubVec1000Inc20(b *testing.B)   { subVecBench(b, 1000, 20) }
 func BenchmarkSubVec10000Inc20(b *testing.B)  { subVecBench(b, 10000, 20) }
 func BenchmarkSubVec100000Inc20(b *testing.B) { subVecBench(b, 100000, 20) }
 func subVecBench(b *testing.B, size, inc int) {
-	x := randVecDense(size, inc, 1, rand.NormFloat64)
-	y := randVecDense(size, inc, 1, rand.NormFloat64)
+	src := rand.NewSource(1)
+	x := randVecDense(size, inc, 1, src)
+	y := randVecDense(size, inc, 1, src)
 	b.ResetTimer()
 	var v VecDense
 	for i := 0; i < b.N; i++ {
@@ -581,14 +585,15 @@ func subVecBench(b *testing.B, size, inc int) {
 	}
 }
 
-func randVecDense(size, inc int, rho float64, rnd func() float64) *VecDense {
+func randVecDense(size, inc int, rho float64, src rand.Source) *VecDense {
 	if size <= 0 {
 		panic("bad vector size")
 	}
+	rnd := rand.New(src)
 	data := make([]float64, size*inc)
 	for i := range data {
-		if rand.Float64() < rho {
-			data[i] = rnd()
+		if rnd.Float64() < rho {
+			data[i] = rnd.NormFloat64()
 		}
 	}
 	return &VecDense{
@@ -605,7 +610,8 @@ func BenchmarkVectorSum100000(b *testing.B) { vectorSumBench(b, 100000) }
 var vectorSumForBench float64
 
 func vectorSumBench(b *testing.B, size int) {
-	a := randVecDense(size, 1, 1.0, rand.NormFloat64)
+	src := rand.NewSource(1)
+	a := randVecDense(size, 1, 1.0, src)
 	b.ResetTimer()
 	for i := 0; i < b.N; i++ {
 		vectorSumForBench = Sum(a)