diff --git a/mat/band_test.go b/mat/band_test.go
index ce6bebfa..1cca02b9 100644
--- a/mat/band_test.go
+++ b/mat/band_test.go
@@ -12,6 +12,7 @@ import (
 )
 
 func TestNewBand(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		data   []float64
 		r, c   int
@@ -159,6 +160,7 @@ func TestNewBand(t *testing.T) {
 }
 
 func TestNewDiagonalRect(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		data  []float64
 		r, c  int
@@ -249,6 +251,7 @@ func TestNewDiagonalRect(t *testing.T) {
 }
 
 func TestBandDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*BandDense{
@@ -302,6 +305,7 @@ func TestBandDenseZero(t *testing.T) {
 }
 
 func TestBandDiagView(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []*BandDense{
 		NewBandDense(1, 1, 0, 0, []float64{1}),
 		NewBandDense(6, 6, 1, 2, []float64{
@@ -326,6 +330,7 @@ func TestBandDiagView(t *testing.T) {
 }
 
 func TestBandAtSet(t *testing.T) {
+	t.Parallel()
 	// 2  3  4  0  0  0
 	// 5  6  7  8  0  0
 	// 0  9 10 11 12  0
diff --git a/mat/cdense_test.go b/mat/cdense_test.go
index 53585064..79862f39 100644
--- a/mat/cdense_test.go
+++ b/mat/cdense_test.go
@@ -7,6 +7,7 @@ package mat
 import "testing"
 
 func TestCDenseNewAtSet(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []struct {
 		a          []complex128
 		rows, cols int
diff --git a/mat/cholesky_test.go b/mat/cholesky_test.go
index 765eb6df..a8325045 100644
--- a/mat/cholesky_test.go
+++ b/mat/cholesky_test.go
@@ -15,6 +15,7 @@ import (
 )
 
 func TestCholesky(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		a *SymDense
 
@@ -74,6 +75,7 @@ func TestCholesky(t *testing.T) {
 }
 
 func TestCholeskyAt(t *testing.T) {
+	t.Parallel()
 	for _, test := range []*SymDense{
 		NewSymDense(3, []float64{
 			53, 59, 37,
@@ -104,6 +106,7 @@ func TestCholeskyAt(t *testing.T) {
 }
 
 func TestCholeskySolveTo(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		a   *SymDense
 		b   *Dense
@@ -151,6 +154,7 @@ func TestCholeskySolveTo(t *testing.T) {
 }
 
 func TestCholeskySolveCholTo(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		a, b *SymDense
 	}{
@@ -218,6 +222,7 @@ func TestCholeskySolveCholTo(t *testing.T) {
 }
 
 func TestCholeskySolveVecTo(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		a   *SymDense
 		b   *VecDense
@@ -265,6 +270,7 @@ func TestCholeskySolveVecTo(t *testing.T) {
 }
 
 func TestCholeskyToSym(t *testing.T) {
+	t.Parallel()
 	for _, test := range []*SymDense{
 		NewSymDense(3, []float64{
 			53, 59, 37,
@@ -287,6 +293,7 @@ func TestCholeskyToSym(t *testing.T) {
 }
 
 func TestCloneCholesky(t *testing.T) {
+	t.Parallel()
 	for _, test := range []*SymDense{
 		NewSymDense(3, []float64{
 			53, 59, 37,
@@ -323,6 +330,7 @@ func TestCloneCholesky(t *testing.T) {
 }
 
 func TestCholeskyInverseTo(t *testing.T) {
+	t.Parallel()
 	for _, n := range []int{1, 3, 5, 9} {
 		data := make([]float64, n*n)
 		for i := range data {
@@ -354,6 +362,7 @@ func TestCholeskyInverseTo(t *testing.T) {
 }
 
 func TestCholeskySymRankOne(t *testing.T) {
+	t.Parallel()
 	rand.Seed(1)
 	for _, n := range []int{1, 2, 3, 4, 5, 7, 10, 20, 50, 100} {
 		for k := 0; k < 50; k++ {
@@ -493,6 +502,7 @@ func TestCholeskySymRankOne(t *testing.T) {
 }
 
 func TestCholeskyExtendVecSym(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []struct {
 		a *SymDense
 	}{
@@ -555,6 +565,7 @@ func TestCholeskyExtendVecSym(t *testing.T) {
 }
 
 func TestCholeskyScale(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []struct {
 		a *SymDense
 		f float64
diff --git a/mat/dense_test.go b/mat/dense_test.go
index f5dbf9f5..46a3c11f 100644
--- a/mat/dense_test.go
+++ b/mat/dense_test.go
@@ -26,6 +26,7 @@ func asBasicTriBanded(t *TriBandDense) TriBanded { return (*basicTriBanded)(t) }
 func asBasicDiagonal(d *DiagDense) Diagonal      { return (*basicDiagonal)(d) }
 
 func TestNewDense(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a          []float64
 		rows, cols int
@@ -168,6 +169,7 @@ func TestNewDense(t *testing.T) {
 }
 
 func TestDenseAtSet(t *testing.T) {
+	t.Parallel()
 	for test, af := range [][][]float64{
 		{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, // even
 		{{1, 2}, {4, 5}, {7, 8}},          // wide
@@ -221,6 +223,7 @@ func TestDenseAtSet(t *testing.T) {
 }
 
 func TestDenseSetRowColumn(t *testing.T) {
+	t.Parallel()
 	for _, as := range [][][]float64{
 		{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
 		{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}},
@@ -259,6 +262,7 @@ func TestDenseSetRowColumn(t *testing.T) {
 }
 
 func TestDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*Dense{
@@ -291,6 +295,7 @@ func TestDenseZero(t *testing.T) {
 }
 
 func TestDenseRowColView(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		mat [][]float64
 	}{
@@ -394,6 +399,7 @@ func TestDenseRowColView(t *testing.T) {
 }
 
 func TestDenseDiagView(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []*Dense{
 		NewDense(1, 1, []float64{1}),
 		NewDense(2, 2, []float64{1, 2, 3, 4}),
@@ -414,6 +420,7 @@ func TestDenseDiagView(t *testing.T) {
 }
 
 func TestDenseGrow(t *testing.T) {
+	t.Parallel()
 	m := &Dense{}
 	m = m.Grow(10, 10).(*Dense)
 	rows, cols := m.Dims()
@@ -469,6 +476,7 @@ func TestDenseGrow(t *testing.T) {
 }
 
 func TestDenseAdd(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b, r [][]float64
 	}{
@@ -556,6 +564,7 @@ func TestDenseAdd(t *testing.T) {
 }
 
 func TestDenseSub(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b, r [][]float64
 	}{
@@ -643,6 +652,7 @@ func TestDenseSub(t *testing.T) {
 }
 
 func TestDenseMulElem(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b, r [][]float64
 	}{
@@ -746,6 +756,7 @@ func (m *Dense) same(b Matrix) bool {
 }
 
 func TestDenseDivElem(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b, r [][]float64
 	}{
@@ -833,6 +844,7 @@ func TestDenseDivElem(t *testing.T) {
 }
 
 func TestDenseMul(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b, r [][]float64
 	}{
@@ -947,6 +959,7 @@ func randDense(size int, rho float64, rnd func() float64) (*Dense, error) {
 }
 
 func TestDenseExp(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a    [][]float64
 		want [][]float64
@@ -1001,6 +1014,7 @@ func TestDenseExp(t *testing.T) {
 }
 
 func TestDensePow(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a    [][]float64
 		n    int
@@ -1088,6 +1102,7 @@ func TestDensePow(t *testing.T) {
 }
 
 func TestDenseKronecker(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b Matrix
 		want *Dense
@@ -1188,6 +1203,7 @@ func TestDenseKronecker(t *testing.T) {
 }
 
 func TestDenseScale(t *testing.T) {
+	t.Parallel()
 	for _, f := range []float64{0.5, 1, 3} {
 		method := func(receiver, a Matrix) {
 			type Scaler interface {
@@ -1204,6 +1220,7 @@ func TestDenseScale(t *testing.T) {
 }
 
 func TestDensePowN(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a   [][]float64
 		mod func(*Dense)
@@ -1244,6 +1261,7 @@ func (m *Dense) iterativePow(a Matrix, n int) {
 }
 
 func TestDenseCloneT(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, want [][]float64
 	}{
@@ -1291,6 +1309,7 @@ func TestDenseCloneT(t *testing.T) {
 }
 
 func TestDenseCopyT(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, want [][]float64
 	}{
@@ -1340,6 +1359,7 @@ func TestDenseCopyT(t *testing.T) {
 }
 
 func TestDenseCopyDenseAlias(t *testing.T) {
+	t.Parallel()
 	for _, trans := range []bool{false, true} {
 		for di := 0; di < 2; di++ {
 			for dj := 0; dj < 2; dj++ {
@@ -1378,6 +1398,7 @@ func TestDenseCopyDenseAlias(t *testing.T) {
 }
 
 func TestDenseCopyVecDenseAlias(t *testing.T) {
+	t.Parallel()
 	for _, horiz := range []bool{false, true} {
 		for do := 0; do < 2; do++ {
 			for di := 0; di < 3; di++ {
@@ -1420,6 +1441,7 @@ func TestDenseCopyVecDenseAlias(t *testing.T) {
 func identity(r, c int, v float64) float64 { return v }
 
 func TestDenseApply(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, want [][]float64
 		fn      func(r, c int, v float64) float64
@@ -1520,6 +1542,7 @@ func TestDenseApply(t *testing.T) {
 }
 
 func TestDenseClone(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a    [][]float64
 		i, j int
@@ -1550,6 +1573,7 @@ func TestDenseClone(t *testing.T) {
 
 // TODO(kortschak) Roll this into testOneInput when it exists.
 func TestDenseCopyPanic(t *testing.T) {
+	t.Parallel()
 	for _, a := range []*Dense{
 		{},
 		{mat: blas64.General{Rows: 1}},
@@ -1571,6 +1595,7 @@ func TestDenseCopyPanic(t *testing.T) {
 }
 
 func TestDenseStack(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b, e [][]float64
 	}{
@@ -1615,6 +1640,7 @@ func TestDenseStack(t *testing.T) {
 }
 
 func TestDenseAugment(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b, e [][]float64
 	}{
@@ -1659,6 +1685,7 @@ func TestDenseAugment(t *testing.T) {
 }
 
 func TestDenseRankOne(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		x     []float64
 		y     []float64
@@ -1737,6 +1764,7 @@ func TestDenseRankOne(t *testing.T) {
 }
 
 func TestDenseOuter(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		x []float64
 		y []float64
@@ -1803,6 +1831,7 @@ func TestDenseOuter(t *testing.T) {
 }
 
 func TestDenseInverse(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a    Matrix
 		want Matrix // nil indicates that a is singular.
@@ -2177,12 +2206,14 @@ func denseMulTransBench(b *testing.B, size int, rho float64) {
 	}
 }
 
-func BenchmarkDenseMulTransDenseSym100Half(b *testing.B)        { denseMulTransSymBench(b, 100, 0.5) }
-func BenchmarkDenseMulTransDenseSym100Tenth(b *testing.B)       { denseMulTransSymBench(b, 100, 0.1) }
-func BenchmarkDenseMulTransDenseSym1000Half(b *testing.B)       { denseMulTransSymBench(b, 1000, 0.5) }
-func BenchmarkDenseMulTransDenseSym1000Tenth(b *testing.B)      { denseMulTransSymBench(b, 1000, 0.1) }
-func BenchmarkDenseMulTransDenseSym1000Hundredth(b *testing.B)  { denseMulTransSymBench(b, 1000, 0.01) }
-func BenchmarkDenseMulTransDenseSym1000Thousandth(b *testing.B) { denseMulTransSymBench(b, 1000, 0.001) }
+func BenchmarkDenseMulTransDenseSym100Half(b *testing.B)       { denseMulTransSymBench(b, 100, 0.5) }
+func BenchmarkDenseMulTransDenseSym100Tenth(b *testing.B)      { denseMulTransSymBench(b, 100, 0.1) }
+func BenchmarkDenseMulTransDenseSym1000Half(b *testing.B)      { denseMulTransSymBench(b, 1000, 0.5) }
+func BenchmarkDenseMulTransDenseSym1000Tenth(b *testing.B)     { denseMulTransSymBench(b, 1000, 0.1) }
+func BenchmarkDenseMulTransDenseSym1000Hundredth(b *testing.B) { denseMulTransSymBench(b, 1000, 0.01) }
+func BenchmarkDenseMulTransDenseSym1000Thousandth(b *testing.B) {
+	denseMulTransSymBench(b, 1000, 0.001)
+}
 func denseMulTransSymBench(b *testing.B, size int, rho float64) {
 	b.StopTimer()
 	a, _ := randDense(size, rho, rand.NormFloat64)
diff --git a/mat/diagonal_test.go b/mat/diagonal_test.go
index e7512df5..9c142b81 100644
--- a/mat/diagonal_test.go
+++ b/mat/diagonal_test.go
@@ -15,6 +15,7 @@ import (
 )
 
 func TestNewDiagDense(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		data  []float64
 		n     int
@@ -59,6 +60,7 @@ func TestNewDiagDense(t *testing.T) {
 }
 
 func TestDiagDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*DiagDense{
@@ -91,6 +93,7 @@ func TestDiagDenseZero(t *testing.T) {
 }
 
 func TestDiagonalStride(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		diag  *DiagDense
 		dense *Dense
@@ -157,6 +160,7 @@ func TestDiagonalStride(t *testing.T) {
 }
 
 func TestDiagFrom(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		mat  Matrix
 		want *Dense
@@ -425,6 +429,7 @@ func testDiagView(t *testing.T, cas int, test diagDenseViewer) {
 }
 
 func TestDiagonalAtSet(t *testing.T) {
+	t.Parallel()
 	for _, n := range []int{1, 3, 8} {
 		for _, nilstart := range []bool{true, false} {
 			var diag *DiagDense
diff --git a/mat/eigen_test.go b/mat/eigen_test.go
index db203acd..769da504 100644
--- a/mat/eigen_test.go
+++ b/mat/eigen_test.go
@@ -14,6 +14,7 @@ import (
 )
 
 func TestEigen(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a *Dense
 
@@ -131,6 +132,7 @@ func cmplxEqualTol(v1, v2 []complex128, tol float64) bool {
 }
 
 func TestSymEigen(t *testing.T) {
+	t.Parallel()
 	// Hand coded tests with results from lapack.
 	for _, test := range []struct {
 		mat *SymDense
diff --git a/mat/errors_test.go b/mat/errors_test.go
index 72b4a5bb..42f93036 100644
--- a/mat/errors_test.go
+++ b/mat/errors_test.go
@@ -16,6 +16,7 @@ func leaksPanic(fn func()) (panicked bool) {
 }
 
 func TestMaybe(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		fn     func()
 		panics bool
diff --git a/mat/format_test.go b/mat/format_test.go
index 444f9219..4efe6629 100644
--- a/mat/format_test.go
+++ b/mat/format_test.go
@@ -11,6 +11,7 @@ import (
 )
 
 func TestFormat(t *testing.T) {
+	t.Parallel()
 	type rp struct {
 		format string
 		output string
diff --git a/mat/gsvd_test.go b/mat/gsvd_test.go
index 2a7afe75..c3dd0348 100644
--- a/mat/gsvd_test.go
+++ b/mat/gsvd_test.go
@@ -5,6 +5,7 @@
 package mat
 
 import (
+	"fmt"
 	"testing"
 
 	"golang.org/x/exp/rand"
@@ -13,8 +14,9 @@ import (
 )
 
 func TestGSVD(t *testing.T) {
+	t.Parallel()
+
 	const tol = 1e-10
-	rnd := rand.New(rand.NewSource(1))
 	for _, test := range []struct {
 		m, p, n int
 	}{
@@ -34,74 +36,80 @@ func TestGSVD(t *testing.T) {
 		m := test.m
 		p := test.p
 		n := test.n
-		for trial := 0; trial < 10; trial++ {
-			a := NewDense(m, n, nil)
-			for i := range a.mat.Data {
-				a.mat.Data[i] = rnd.NormFloat64()
-			}
-			aCopy := DenseCopyOf(a)
+		t.Run(fmt.Sprintf("%v", test), func(t *testing.T) {
+			t.Parallel()
 
-			b := NewDense(p, n, nil)
-			for i := range b.mat.Data {
-				b.mat.Data[i] = rnd.NormFloat64()
-			}
-			bCopy := DenseCopyOf(b)
+			rnd := rand.New(rand.NewSource(1))
+			for trial := 0; trial < 10; trial++ {
+				a := NewDense(m, n, nil)
+				for i := range a.mat.Data {
+					a.mat.Data[i] = rnd.NormFloat64()
+				}
+				aCopy := DenseCopyOf(a)
 
-			// Test Full decomposition.
-			var gsvd GSVD
-			ok := gsvd.Factorize(a, b, GSVDU|GSVDV|GSVDQ)
-			if !ok {
-				t.Errorf("GSVD factorization failed")
-			}
-			if !Equal(a, aCopy) {
-				t.Errorf("A changed during call to GSVD.Factorize with GSVDU|GSVDV|GSVDQ")
-			}
-			if !Equal(b, bCopy) {
-				t.Errorf("B changed during call to GSVD.Factorize with GSVDU|GSVDV|GSVDQ")
-			}
-			c, s, sigma1, sigma2, zeroR, u, v, q := extractGSVD(&gsvd)
-			var ansU, ansV, d1R, d2R Dense
-			ansU.Product(u.T(), a, q)
-			ansV.Product(v.T(), b, q)
-			d1R.Mul(sigma1, zeroR)
-			d2R.Mul(sigma2, zeroR)
-			if !EqualApprox(&ansU, &d1R, tol) {
-				t.Errorf("Answer mismatch with GSVDU|GSVDV|GSVDQ\nUᵀ * A * Q:\n% 0.2f\nΣ₁ * [ 0 R ]:\n% 0.2f",
-					Formatted(&ansU), Formatted(&d1R))
-			}
-			if !EqualApprox(&ansV, &d2R, tol) {
-				t.Errorf("Answer mismatch with GSVDU|GSVDV|GSVDQ\nVᵀ * B  *Q:\n% 0.2f\nΣ₂ * [ 0 R ]:\n% 0.2f",
-					Formatted(&d2R), Formatted(&ansV))
-			}
+				b := NewDense(p, n, nil)
+				for i := range b.mat.Data {
+					b.mat.Data[i] = rnd.NormFloat64()
+				}
+				bCopy := DenseCopyOf(b)
 
-			// Check C^2 + S^2 = I.
-			for i := range c {
-				d := c[i]*c[i] + s[i]*s[i]
-				if !floats.EqualWithinAbsOrRel(d, 1, 1e-14, 1e-14) {
-					t.Errorf("c_%d^2 + s_%d^2 != 1: got: %v", i, i, d)
+				// Test Full decomposition.
+				var gsvd GSVD
+				ok := gsvd.Factorize(a, b, GSVDU|GSVDV|GSVDQ)
+				if !ok {
+					t.Errorf("GSVD factorization failed")
+				}
+				if !Equal(a, aCopy) {
+					t.Errorf("A changed during call to GSVD.Factorize with GSVDU|GSVDV|GSVDQ")
+				}
+				if !Equal(b, bCopy) {
+					t.Errorf("B changed during call to GSVD.Factorize with GSVDU|GSVDV|GSVDQ")
+				}
+				c, s, sigma1, sigma2, zeroR, u, v, q := extractGSVD(&gsvd)
+				var ansU, ansV, d1R, d2R Dense
+				ansU.Product(u.T(), a, q)
+				ansV.Product(v.T(), b, q)
+				d1R.Mul(sigma1, zeroR)
+				d2R.Mul(sigma2, zeroR)
+				if !EqualApprox(&ansU, &d1R, tol) {
+					t.Errorf("Answer mismatch with GSVDU|GSVDV|GSVDQ\nUᵀ * A * Q:\n% 0.2f\nΣ₁ * [ 0 R ]:\n% 0.2f",
+						Formatted(&ansU), Formatted(&d1R))
+				}
+				if !EqualApprox(&ansV, &d2R, tol) {
+					t.Errorf("Answer mismatch with GSVDU|GSVDV|GSVDQ\nVᵀ * B  *Q:\n% 0.2f\nΣ₂ * [ 0 R ]:\n% 0.2f",
+						Formatted(&d2R), Formatted(&ansV))
+				}
+
+				// Check C^2 + S^2 = I.
+				for i := range c {
+					d := c[i]*c[i] + s[i]*s[i]
+					if !floats.EqualWithinAbsOrRel(d, 1, 1e-14, 1e-14) {
+						t.Errorf("c_%d^2 + s_%d^2 != 1: got: %v", i, i, d)
+					}
+				}
+
+				// Test None decomposition.
+				ok = gsvd.Factorize(a, b, GSVDNone)
+				if !ok {
+					t.Errorf("GSVD factorization failed")
+				}
+				if !Equal(a, aCopy) {
+					t.Errorf("A changed during call to GSVD with GSVDNone")
+				}
+				if !Equal(b, bCopy) {
+					t.Errorf("B changed during call to GSVD with GSVDNone")
+				}
+				cNone := gsvd.ValuesA(nil)
+				if !floats.EqualApprox(c, cNone, tol) {
+					t.Errorf("Singular value mismatch between GSVDU|GSVDV|GSVDQ and GSVDNone decomposition")
+				}
+				sNone := gsvd.ValuesB(nil)
+				if !floats.EqualApprox(s, sNone, tol) {
+					t.Errorf("Singular value mismatch between GSVDU|GSVDV|GSVDQ and GSVDNone decomposition")
 				}
 			}
+		})
 
-			// Test None decomposition.
-			ok = gsvd.Factorize(a, b, GSVDNone)
-			if !ok {
-				t.Errorf("GSVD factorization failed")
-			}
-			if !Equal(a, aCopy) {
-				t.Errorf("A changed during call to GSVD with GSVDNone")
-			}
-			if !Equal(b, bCopy) {
-				t.Errorf("B changed during call to GSVD with GSVDNone")
-			}
-			cNone := gsvd.ValuesA(nil)
-			if !floats.EqualApprox(c, cNone, tol) {
-				t.Errorf("Singular value mismatch between GSVDU|GSVDV|GSVDQ and GSVDNone decomposition")
-			}
-			sNone := gsvd.ValuesB(nil)
-			if !floats.EqualApprox(s, sNone, tol) {
-				t.Errorf("Singular value mismatch between GSVDU|GSVDV|GSVDQ and GSVDNone decomposition")
-			}
-		}
 	}
 }
 
diff --git a/mat/hogsvd_test.go b/mat/hogsvd_test.go
index 406cef1c..07886ddd 100644
--- a/mat/hogsvd_test.go
+++ b/mat/hogsvd_test.go
@@ -11,6 +11,7 @@ import (
 )
 
 func TestHOGSVD(t *testing.T) {
+	t.Parallel()
 	const tol = 1e-10
 	rnd := rand.New(rand.NewSource(1))
 	for cas, test := range []struct {
diff --git a/mat/inner_test.go b/mat/inner_test.go
index 3db9d28f..ad53076e 100644
--- a/mat/inner_test.go
+++ b/mat/inner_test.go
@@ -13,6 +13,7 @@ import (
 )
 
 func TestInner(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		x []float64
 		y []float64
@@ -92,6 +93,7 @@ func TestInner(t *testing.T) {
 }
 
 func TestInnerSym(t *testing.T) {
+	t.Parallel()
 	for _, inc := range []struct{ x, y int }{
 		{1, 1},
 		{1, 2},
diff --git a/mat/io_test.go b/mat/io_test.go
index 28b26aec..120d5925 100644
--- a/mat/io_test.go
+++ b/mat/io_test.go
@@ -92,6 +92,7 @@ var denseData = []struct {
 }
 
 func TestDenseMarshal(t *testing.T) {
+	t.Parallel()
 	for i, test := range denseData {
 		buf, err := test.want.MarshalBinary()
 		if err != nil {
@@ -117,6 +118,7 @@ func TestDenseMarshal(t *testing.T) {
 }
 
 func TestDenseMarshalTo(t *testing.T) {
+	t.Parallel()
 	for i, test := range denseData {
 		buf := new(bytes.Buffer)
 		n, err := test.want.MarshalBinaryTo(buf)
@@ -143,6 +145,7 @@ func TestDenseMarshalTo(t *testing.T) {
 }
 
 func TestDenseUnmarshal(t *testing.T) {
+	t.Parallel()
 	for i, test := range denseData {
 		var v Dense
 		err := v.UnmarshalBinary(test.raw)
@@ -163,6 +166,7 @@ func TestDenseUnmarshal(t *testing.T) {
 }
 
 func TestDenseUnmarshalFrom(t *testing.T) {
+	t.Parallel()
 	for i, test := range denseData {
 		var v Dense
 		buf := bytes.NewReader(test.raw)
@@ -189,6 +193,7 @@ func TestDenseUnmarshalFrom(t *testing.T) {
 }
 
 func TestDenseUnmarshalFromError(t *testing.T) {
+	t.Parallel()
 	test := denseData[1]
 	for i, tt := range []struct {
 		beg int
@@ -257,6 +262,7 @@ func TestDenseUnmarshalFromError(t *testing.T) {
 }
 
 func TestDenseIORoundTrip(t *testing.T) {
+	t.Parallel()
 	for i, test := range denseData {
 		buf, err := test.want.MarshalBinary()
 		if err != nil {
@@ -373,6 +379,7 @@ var vectorData = []struct {
 }
 
 func TestVecDenseMarshal(t *testing.T) {
+	t.Parallel()
 	for i, test := range vectorData {
 		buf, err := test.want.MarshalBinary()
 		if err != nil {
@@ -398,6 +405,7 @@ func TestVecDenseMarshal(t *testing.T) {
 }
 
 func TestVecDenseMarshalTo(t *testing.T) {
+	t.Parallel()
 	for i, test := range vectorData {
 		buf := new(bytes.Buffer)
 		n, err := test.want.MarshalBinaryTo(buf)
@@ -424,6 +432,7 @@ func TestVecDenseMarshalTo(t *testing.T) {
 }
 
 func TestVecDenseUnmarshal(t *testing.T) {
+	t.Parallel()
 	for i, test := range vectorData {
 		var v VecDense
 		err := v.UnmarshalBinary(test.raw)
@@ -444,6 +453,7 @@ func TestVecDenseUnmarshal(t *testing.T) {
 }
 
 func TestVecDenseUnmarshalFrom(t *testing.T) {
+	t.Parallel()
 	for i, test := range vectorData {
 		var v VecDense
 		buf := bytes.NewReader(test.raw)
@@ -472,6 +482,7 @@ func TestVecDenseUnmarshalFrom(t *testing.T) {
 }
 
 func TestVecDenseUnmarshalFromError(t *testing.T) {
+	t.Parallel()
 	test := vectorData[1]
 	for i, tt := range []struct {
 		beg int
@@ -528,6 +539,7 @@ func TestVecDenseUnmarshalFromError(t *testing.T) {
 }
 
 func TestVecDenseIORoundTrip(t *testing.T) {
+	t.Parallel()
 	for i, test := range vectorData {
 		buf, err := test.want.MarshalBinary()
 		if err != nil {
diff --git a/mat/lq_test.go b/mat/lq_test.go
index f76be059..3fcc6ef0 100644
--- a/mat/lq_test.go
+++ b/mat/lq_test.go
@@ -11,6 +11,7 @@ import (
 )
 
 func TestLQ(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		m, n int
 	}{
@@ -48,6 +49,7 @@ func TestLQ(t *testing.T) {
 }
 
 func TestLQSolveTo(t *testing.T) {
+	t.Parallel()
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n, bc int
@@ -109,6 +111,7 @@ func TestLQSolveTo(t *testing.T) {
 }
 
 func TestLQSolveToVec(t *testing.T) {
+	t.Parallel()
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n int
@@ -165,6 +168,7 @@ func TestLQSolveToVec(t *testing.T) {
 }
 
 func TestLQSolveToCond(t *testing.T) {
+	t.Parallel()
 	for _, test := range []*Dense{
 		NewDense(2, 2, []float64{1, 0, 0, 1e-20}),
 		NewDense(2, 3, []float64{1, 0, 0, 0, 1e-20, 0}),
diff --git a/mat/lu_test.go b/mat/lu_test.go
index ec2d25c4..602b4056 100644
--- a/mat/lu_test.go
+++ b/mat/lu_test.go
@@ -11,6 +11,7 @@ import (
 )
 
 func TestLUD(t *testing.T) {
+	t.Parallel()
 	for _, n := range []int{1, 5, 10, 11, 50} {
 		a := NewDense(n, n, nil)
 		for i := 0; i < n; i++ {
@@ -39,6 +40,7 @@ func TestLUD(t *testing.T) {
 }
 
 func TestLURankOne(t *testing.T) {
+	t.Parallel()
 	for _, pivoting := range []bool{true} {
 		for _, n := range []int{3, 10, 50} {
 			// Construct a random LU factorization
@@ -106,6 +108,7 @@ func luReconstruct(lu *LU) *Dense {
 }
 
 func TestLUSolveTo(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		n, bc int
 	}{
@@ -143,6 +146,7 @@ func TestLUSolveTo(t *testing.T) {
 }
 
 func TestLUSolveToCond(t *testing.T) {
+	t.Parallel()
 	for _, test := range []*Dense{
 		NewDense(2, 2, []float64{1, 0, 0, 1e-20}),
 	} {
@@ -164,6 +168,7 @@ func TestLUSolveToCond(t *testing.T) {
 }
 
 func TestLUSolveVecTo(t *testing.T) {
+	t.Parallel()
 	for _, n := range []int{5, 10} {
 		a := NewDense(n, n, nil)
 		for i := 0; i < n; i++ {
diff --git a/mat/matrix_test.go b/mat/matrix_test.go
index 1d9ae570..ab75df5f 100644
--- a/mat/matrix_test.go
+++ b/mat/matrix_test.go
@@ -60,6 +60,7 @@ func eye(n int) *Dense {
 }
 
 func TestCol(t *testing.T) {
+	t.Parallel()
 	for id, af := range [][][]float64{
 		{
 			{1, 2, 3},
@@ -133,6 +134,7 @@ func TestCol(t *testing.T) {
 }
 
 func TestRow(t *testing.T) {
+	t.Parallel()
 	for id, af := range [][][]float64{
 		{
 			{1, 2, 3},
@@ -201,6 +203,7 @@ func TestRow(t *testing.T) {
 }
 
 func TestCond(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a       *Dense
 		condOne float64
@@ -290,6 +293,7 @@ func TestCond(t *testing.T) {
 }
 
 func TestDet(t *testing.T) {
+	t.Parallel()
 	for c, test := range []struct {
 		a   *Dense
 		ans float64
@@ -411,6 +415,7 @@ func (v *rawVector) RawVector() blas64.Vector {
 }
 
 func TestDot(t *testing.T) {
+	t.Parallel()
 	f := func(a, b Matrix) interface{} {
 		return Dot(a.(Vector), b.(Vector))
 	}
@@ -432,6 +437,7 @@ func TestDot(t *testing.T) {
 }
 
 func TestEqual(t *testing.T) {
+	t.Parallel()
 	f := func(a, b Matrix) interface{} {
 		return Equal(a, b)
 	}
@@ -442,6 +448,7 @@ func TestEqual(t *testing.T) {
 }
 
 func TestMax(t *testing.T) {
+	t.Parallel()
 	// A direct test of Max with *Dense arguments is in TestNewDense.
 	f := func(a Matrix) interface{} {
 		return Max(a)
@@ -453,6 +460,7 @@ func TestMax(t *testing.T) {
 }
 
 func TestMin(t *testing.T) {
+	t.Parallel()
 	// A direct test of Min with *Dense arguments is in TestNewDense.
 	f := func(a Matrix) interface{} {
 		return Min(a)
@@ -464,6 +472,7 @@ func TestMin(t *testing.T) {
 }
 
 func TestNorm(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a    [][]float64
 		ord  float64
@@ -520,6 +529,7 @@ func TestNorm(t *testing.T) {
 }
 
 func TestNormZero(t *testing.T) {
+	t.Parallel()
 	for _, a := range []Matrix{
 		&Dense{},
 		&SymDense{},
@@ -542,6 +552,7 @@ func TestNormZero(t *testing.T) {
 }
 
 func TestSum(t *testing.T) {
+	t.Parallel()
 	f := func(a Matrix) interface{} {
 		return Sum(a)
 	}
@@ -552,6 +563,7 @@ func TestSum(t *testing.T) {
 }
 
 func TestTrace(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		a     *Dense
 		trace float64
@@ -576,6 +588,7 @@ func TestTrace(t *testing.T) {
 }
 
 func TestTracer(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		a    Tracer
 		want float64
@@ -617,6 +630,7 @@ func TestTracer(t *testing.T) {
 }
 
 func TestDoer(t *testing.T) {
+	t.Parallel()
 	type MatrixDoer interface {
 		Matrix
 		NonZeroDoer
@@ -708,6 +722,7 @@ func TestDoer(t *testing.T) {
 }
 
 func TestMulVecToer(t *testing.T) {
+	t.Parallel()
 	const tol = 1e-14
 
 	rnd := rand.New(rand.NewSource(1))
diff --git a/mat/mul_test.go b/mat/mul_test.go
index d279aed3..5dafeac1 100644
--- a/mat/mul_test.go
+++ b/mat/mul_test.go
@@ -16,6 +16,7 @@ import (
 
 // TODO: Need to add tests where one is overwritten.
 func TestMulTypes(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		ar     int
 		ac     int
diff --git a/mat/pool_test.go b/mat/pool_test.go
index e66142d7..196d24ef 100644
--- a/mat/pool_test.go
+++ b/mat/pool_test.go
@@ -13,6 +13,7 @@ import (
 )
 
 func TestPool(t *testing.T) {
+	t.Parallel()
 	for i := 1; i < 10; i++ {
 		for j := 1; j < 10; j++ {
 			m := NewDense(i, j, nil)
diff --git a/mat/product_test.go b/mat/product_test.go
index 31b18d1b..bf659099 100644
--- a/mat/product_test.go
+++ b/mat/product_test.go
@@ -91,6 +91,7 @@ var productTests = []struct {
 }
 
 func TestProduct(t *testing.T) {
+	t.Parallel()
 	for _, test := range productTests {
 		dimensions := test.factors
 		if dimensions == nil && test.n > 0 {
diff --git a/mat/qr_test.go b/mat/qr_test.go
index b800872c..4cbd173d 100644
--- a/mat/qr_test.go
+++ b/mat/qr_test.go
@@ -14,6 +14,7 @@ import (
 )
 
 func TestQR(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		m, n int
 	}{
@@ -72,6 +73,7 @@ func isOrthonormal(q *Dense, tol float64) bool {
 }
 
 func TestQRSolveTo(t *testing.T) {
+	t.Parallel()
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n, bc int
@@ -133,6 +135,7 @@ func TestQRSolveTo(t *testing.T) {
 }
 
 func TestQRSolveVecTo(t *testing.T) {
+	t.Parallel()
 	for _, trans := range []bool{false, true} {
 		for _, test := range []struct {
 			m, n int
@@ -189,6 +192,7 @@ func TestQRSolveVecTo(t *testing.T) {
 }
 
 func TestQRSolveCondTo(t *testing.T) {
+	t.Parallel()
 	for _, test := range []*Dense{
 		NewDense(2, 2, []float64{1, 0, 0, 1e-20}),
 		NewDense(3, 2, []float64{1, 0, 0, 1e-20, 0, 0}),
diff --git a/mat/shadow_test.go b/mat/shadow_test.go
index ab0f78d9..8c7b8475 100644
--- a/mat/shadow_test.go
+++ b/mat/shadow_test.go
@@ -11,6 +11,7 @@ import (
 )
 
 func TestDenseOverlaps(t *testing.T) {
+	t.Parallel()
 	type view struct {
 		i, j, r, c int
 		*Dense
@@ -91,6 +92,7 @@ func intervalsOverlap(a, b interval) bool {
 
 // See https://github.com/gonum/matrix/issues/359 for details.
 func TestIssue359(t *testing.T) {
+	t.Parallel()
 	for xi := 0; xi < 2; xi++ {
 		for xj := 0; xj < 2; xj++ {
 			for yi := 0; yi < 2; yi++ {
diff --git a/mat/solve_test.go b/mat/solve_test.go
index 985af8b9..bb28532a 100644
--- a/mat/solve_test.go
+++ b/mat/solve_test.go
@@ -11,6 +11,7 @@ import (
 )
 
 func TestSolve(t *testing.T) {
+	t.Parallel()
 	// Hand-coded cases.
 	for _, test := range []struct {
 		a         [][]float64
@@ -245,6 +246,7 @@ func TestSolve(t *testing.T) {
 }
 
 func TestSolveVec(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		m, n int
 	}{
diff --git a/mat/svd_test.go b/mat/svd_test.go
index 6349eeca..c5d62dba 100644
--- a/mat/svd_test.go
+++ b/mat/svd_test.go
@@ -13,6 +13,7 @@ import (
 )
 
 func TestSVD(t *testing.T) {
+	t.Parallel()
 	// Hand coded tests
 	for _, test := range []struct {
 		a *Dense
diff --git a/mat/symband_test.go b/mat/symband_test.go
index 6484c6fe..3a757227 100644
--- a/mat/symband_test.go
+++ b/mat/symband_test.go
@@ -13,6 +13,7 @@ import (
 )
 
 func TestNewSymBand(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		data  []float64
 		n     int
@@ -79,6 +80,7 @@ func TestNewSymBand(t *testing.T) {
 }
 
 func TestSymBandAtSet(t *testing.T) {
+	t.Parallel()
 	// 1  2  3  0  0  0
 	// 2  4  5  6  0  0
 	// 3  5  7  8  9  0
@@ -181,6 +183,7 @@ func TestSymBandAtSet(t *testing.T) {
 }
 
 func TestSymBandDiagView(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []*SymBandDense{
 		NewSymBandDense(1, 0, []float64{1}),
 		NewSymBandDense(6, 2, []float64{
@@ -197,6 +200,7 @@ func TestSymBandDiagView(t *testing.T) {
 }
 
 func TestSymBandDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*SymBandDense{
diff --git a/mat/symmetric_test.go b/mat/symmetric_test.go
index f48e9555..867e47e3 100644
--- a/mat/symmetric_test.go
+++ b/mat/symmetric_test.go
@@ -18,6 +18,7 @@ import (
 )
 
 func TestNewSymmetric(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		data []float64
 		n    int
@@ -67,6 +68,7 @@ func TestNewSymmetric(t *testing.T) {
 }
 
 func TestSymAtSet(t *testing.T) {
+	t.Parallel()
 	sym := &SymDense{
 		mat: blas64.Symmetric{
 			N:      3,
@@ -130,6 +132,7 @@ func TestSymAtSet(t *testing.T) {
 }
 
 func TestSymDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*SymDense{
@@ -162,6 +165,7 @@ func TestSymDenseZero(t *testing.T) {
 }
 
 func TestSymDiagView(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []*SymDense{
 		NewSymDense(1, []float64{1}),
 		NewSymDense(2, []float64{1, 2, 2, 3}),
@@ -172,6 +176,7 @@ func TestSymDiagView(t *testing.T) {
 }
 
 func TestSymAdd(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		n int
 	}{
@@ -233,6 +238,7 @@ func TestSymAdd(t *testing.T) {
 }
 
 func TestCopy(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		n int
 	}{
@@ -264,6 +270,7 @@ func TestCopy(t *testing.T) {
 // TODO(kortschak) Roll this into testOneInput when it exists.
 // https://github.com/gonum/matrix/issues/171
 func TestSymCopyPanic(t *testing.T) {
+	t.Parallel()
 	var (
 		a SymDense
 		n int
@@ -279,6 +286,7 @@ func TestSymCopyPanic(t *testing.T) {
 }
 
 func TestSymRankOne(t *testing.T) {
+	t.Parallel()
 	const tol = 1e-15
 
 	for _, test := range []struct {
@@ -365,6 +373,7 @@ func TestSymRankOne(t *testing.T) {
 }
 
 func TestIssue250SymRankOne(t *testing.T) {
+	t.Parallel()
 	x := NewVecDense(5, []float64{1, 2, 3, 4, 5})
 	var s1, s2 SymDense
 	s1.SymRankOne(NewSymDense(5, nil), 1, x)
@@ -376,6 +385,7 @@ func TestIssue250SymRankOne(t *testing.T) {
 }
 
 func TestRankTwo(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		n int
 	}{
@@ -434,6 +444,7 @@ func TestRankTwo(t *testing.T) {
 }
 
 func TestSymRankK(t *testing.T) {
+	t.Parallel()
 	alpha := 3.0
 	method := func(receiver, a, b Matrix) {
 		type SymRankKer interface {
@@ -462,6 +473,7 @@ func TestSymRankK(t *testing.T) {
 }
 
 func TestSymOuterK(t *testing.T) {
+	t.Parallel()
 	for _, f := range []float64{0.5, 1, 3} {
 		method := func(receiver, x Matrix) {
 			type SymOuterKer interface {
@@ -479,6 +491,7 @@ func TestSymOuterK(t *testing.T) {
 }
 
 func TestIssue250SymOuterK(t *testing.T) {
+	t.Parallel()
 	x := NewVecDense(5, []float64{1, 2, 3, 4, 5})
 	var s1, s2 SymDense
 	s1.SymOuterK(1, x)
@@ -490,6 +503,7 @@ func TestIssue250SymOuterK(t *testing.T) {
 }
 
 func TestScaleSym(t *testing.T) {
+	t.Parallel()
 	for _, f := range []float64{0.5, 1, 3} {
 		method := func(receiver, a Matrix) {
 			type ScaleSymer interface {
@@ -506,6 +520,7 @@ func TestScaleSym(t *testing.T) {
 }
 
 func TestSubsetSym(t *testing.T) {
+	t.Parallel()
 	for _, test := range []struct {
 		a    *SymDense
 		dims []int
@@ -587,6 +602,7 @@ func TestSubsetSym(t *testing.T) {
 }
 
 func TestViewGrowSquare(t *testing.T) {
+	t.Parallel()
 	// n is the size of the original SymDense.
 	// The first view uses start1, span1. The second view uses start2, span2 on
 	// the first view.
@@ -669,6 +685,7 @@ func TestViewGrowSquare(t *testing.T) {
 }
 
 func TestPowPSD(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []struct {
 		a   *SymDense
 		pow float64
diff --git a/mat/triangular_test.go b/mat/triangular_test.go
index 38682bc9..191babf0 100644
--- a/mat/triangular_test.go
+++ b/mat/triangular_test.go
@@ -17,6 +17,7 @@ import (
 )
 
 func TestNewTriangular(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		data []float64
 		n    int
@@ -66,6 +67,7 @@ func TestNewTriangular(t *testing.T) {
 }
 
 func TestTriAtSet(t *testing.T) {
+	t.Parallel()
 	tri := &TriDense{
 		mat: blas64.Triangular{
 			N:      3,
@@ -141,6 +143,7 @@ func TestTriAtSet(t *testing.T) {
 }
 
 func TestTriDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*TriDense{
@@ -186,6 +189,7 @@ func TestTriDenseZero(t *testing.T) {
 }
 
 func TestTriDiagView(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []*TriDense{
 		NewTriDense(1, Upper, []float64{1}),
 		NewTriDense(2, Upper, []float64{1, 2, 0, 3}),
@@ -199,6 +203,7 @@ func TestTriDiagView(t *testing.T) {
 }
 
 func TestTriDenseCopy(t *testing.T) {
+	t.Parallel()
 	for i := 0; i < 100; i++ {
 		size := rand.Intn(100)
 		r, err := randDense(size, 0.9, rand.NormFloat64)
@@ -249,6 +254,7 @@ func TestTriDenseCopy(t *testing.T) {
 }
 
 func TestTriTriDenseCopy(t *testing.T) {
+	t.Parallel()
 	for i := 0; i < 100; i++ {
 		size := rand.Intn(100)
 		r, err := randDense(size, 1, rand.NormFloat64)
@@ -301,6 +307,7 @@ func TestTriTriDenseCopy(t *testing.T) {
 }
 
 func TestTriInverse(t *testing.T) {
+	t.Parallel()
 	for _, kind := range []TriKind{Upper, Lower} {
 		for _, n := range []int{1, 3, 5, 9} {
 			data := make([]float64, n*n)
@@ -325,6 +332,7 @@ func TestTriInverse(t *testing.T) {
 }
 
 func TestTriMul(t *testing.T) {
+	t.Parallel()
 	method := func(receiver, a, b Matrix) {
 		type MulTrier interface {
 			MulTri(a, b Triangular)
@@ -380,6 +388,7 @@ func TestTriMul(t *testing.T) {
 }
 
 func TestScaleTri(t *testing.T) {
+	t.Parallel()
 	for _, f := range []float64{0.5, 1, 3} {
 		method := func(receiver, a Matrix) {
 			type ScaleTrier interface {
@@ -397,6 +406,7 @@ func TestScaleTri(t *testing.T) {
 }
 
 func TestCopySymIntoTriangle(t *testing.T) {
+	t.Parallel()
 	nan := math.NaN()
 	for tc, test := range []struct {
 		n     int
@@ -515,6 +525,7 @@ func TestCopySymIntoTriangle(t *testing.T) {
 }
 
 func TestTriSliceTri(t *testing.T) {
+	t.Parallel()
 	rnd := rand.New(rand.NewSource(1))
 	for cas, test := range []struct {
 		n, start1, span1, start2, span2 int
diff --git a/mat/triband_test.go b/mat/triband_test.go
index 4bc97a15..6757a698 100644
--- a/mat/triband_test.go
+++ b/mat/triband_test.go
@@ -13,6 +13,7 @@ import (
 )
 
 func TestNewTriBand(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []struct {
 		data  []float64
 		n, k  int
@@ -215,6 +216,7 @@ func TestNewTriBand(t *testing.T) {
 }
 
 func TestTriBandAtSetUpper(t *testing.T) {
+	t.Parallel()
 	for _, kind := range []TriKind{Upper, Lower} {
 		var band *TriBandDense
 		var data []float64
@@ -334,6 +336,7 @@ func TestTriBandAtSetUpper(t *testing.T) {
 }
 
 func TestTriBandDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*TriBandDense{
@@ -385,6 +388,7 @@ func TestTriBandDenseZero(t *testing.T) {
 }
 
 func TestTriBandDiagView(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []*TriBandDense{
 		NewTriBandDense(1, 0, Upper, []float64{1}),
 		NewTriBandDense(4, 0, Upper, []float64{1, 2, 3, 4}),
diff --git a/mat/vector_test.go b/mat/vector_test.go
index 28a67505..75842806 100644
--- a/mat/vector_test.go
+++ b/mat/vector_test.go
@@ -14,6 +14,7 @@ import (
 )
 
 func TestNewVecDense(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		n      int
 		data   []float64
@@ -57,6 +58,7 @@ func TestNewVecDense(t *testing.T) {
 }
 
 func TestCap(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		vector *VecDense
 		want   int
@@ -121,6 +123,7 @@ func TestCap(t *testing.T) {
 }
 
 func TestVecDenseAtSet(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		vector *VecDense
 	}{
@@ -182,6 +185,7 @@ func TestVecDenseAtSet(t *testing.T) {
 }
 
 func TestVecDenseZero(t *testing.T) {
+	t.Parallel()
 	// Elements that equal 1 should be set to zero, elements that equal -1
 	// should remain unchanged.
 	for _, test := range []*VecDense{
@@ -214,6 +218,7 @@ func TestVecDenseZero(t *testing.T) {
 }
 
 func TestVecDenseMul(t *testing.T) {
+	t.Parallel()
 	method := func(receiver, a, b Matrix) {
 		type mulVecer interface {
 			MulVec(a Matrix, b Vector)
@@ -237,6 +242,7 @@ func TestVecDenseMul(t *testing.T) {
 }
 
 func TestVecDenseScale(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a     Vector
 		alpha float64
@@ -311,6 +317,7 @@ func TestVecDenseScale(t *testing.T) {
 }
 
 func TestCopyVec(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		src   *VecDense
 		dst   *VecDense
@@ -337,6 +344,7 @@ func TestCopyVec(t *testing.T) {
 }
 
 func TestVecDenseAddScaled(t *testing.T) {
+	t.Parallel()
 	for _, alpha := range []float64{0, 1, -1, 2.3, -2.3} {
 		method := func(receiver, a, b Matrix) {
 			type addScaledVecer interface {
@@ -355,6 +363,7 @@ func TestVecDenseAddScaled(t *testing.T) {
 }
 
 func TestVecDenseAdd(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b Vector
 		want *VecDense
@@ -384,6 +393,7 @@ func TestVecDenseAdd(t *testing.T) {
 }
 
 func TestVecDenseSub(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b Vector
 		want *VecDense
@@ -413,6 +423,7 @@ func TestVecDenseSub(t *testing.T) {
 }
 
 func TestVecDenseMulElem(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b Vector
 		want *VecDense
@@ -442,6 +453,7 @@ func TestVecDenseMulElem(t *testing.T) {
 }
 
 func TestVecDenseDivElem(t *testing.T) {
+	t.Parallel()
 	for i, test := range []struct {
 		a, b Vector
 		want *VecDense