diff --git a/mat/band.go b/mat/band.go
index d196a5ae..5728f48d 100644
--- a/mat/band.go
+++ b/mat/band.go
@@ -107,6 +107,7 @@ func (t TransposeBand) UntransposeBand() Banded {
 // BandDense will be reflected in data. If neither of these is true, NewBandDense
 // will panic. kl must be at least zero and less r, and ku must be at least zero and
 // less than c, otherwise NewBandDense will panic.
+// NewBandDense will panic if either r or c is zero.
 //
 // The data must be arranged in row-major order constructed by removing the zeros
 // from the rows outside the band and aligning the diagonals. For example, the matrix
@@ -126,7 +127,10 @@ func (t TransposeBand) UntransposeBand() Banded {
 // which is passed to NewBandDense as []float64{*, 1, 2, 3, 4, ...} with kl=1 and ku=2.
 // Only the values in the band portion of the matrix are used.
 func NewBandDense(r, c, kl, ku int, data []float64) *BandDense {
-	if r < 0 || c < 0 || kl < 0 || ku < 0 {
+	if r <= 0 || c <= 0 || kl < 0 || ku < 0 {
+		if r == 0 || c == 0 {
+			panic(ErrZeroLength)
+		}
 		panic("mat: negative dimension")
 	}
 	if kl+1 > r || ku+1 > c {
diff --git a/mat/dense.go b/mat/dense.go
index 2057ceea..b21652ab 100644
--- a/mat/dense.go
+++ b/mat/dense.go
@@ -38,11 +38,15 @@ type Dense struct {
 // a new slice is allocated for the backing slice. If len(data) == r*c, data is
 // used as the backing slice, and changes to the elements of the returned Dense
 // will be reflected in data. If neither of these is true, NewDense will panic.
+// NewDense will panic if either r or c is zero.
 //
 // The data must be arranged in row-major order, i.e. the (i*c + j)-th
 // element in the data slice is the {i, j}-th element in the matrix.
 func NewDense(r, c int, data []float64) *Dense {
-	if r < 0 || c < 0 {
+	if r <= 0 || c <= 0 {
+		if r == 0 || c == 0 {
+			panic(ErrZeroLength)
+		}
 		panic("mat: negative dimension")
 	}
 	if data != nil && r*c != len(data) {
diff --git a/mat/diagonal.go b/mat/diagonal.go
index 3dbaaf63..acd92d34 100644
--- a/mat/diagonal.go
+++ b/mat/diagonal.go
@@ -42,9 +42,12 @@ type DiagDense struct {
 
 // NewDiagonal creates a new Diagonal matrix with n rows and n columns.
 // The length of data must be n or data must be nil, otherwise NewDiagonal
-// will panic.
+// will panic. NewDiagonal will panic if n is zero.
 func NewDiagonal(n int, data []float64) *DiagDense {
-	if n < 0 {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
 		panic("mat: negative dimension")
 	}
 	if data == nil {
diff --git a/mat/io_test.go b/mat/io_test.go
index 22e18f4b..6589590c 100644
--- a/mat/io_test.go
+++ b/mat/io_test.go
@@ -30,7 +30,7 @@ var denseData = []struct {
 }{
 	{
 		raw:  []byte("\x01\x00\x00\x00GFA\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00"),
-		want: NewDense(0, 0, []float64{}),
+		want: &Dense{},
 		err:  ErrZeroLength,
 		eq:   Equal,
 	},
@@ -305,7 +305,7 @@ var vectorData = []struct {
 }{
 	{
 		raw:  []byte("\x01\x00\x00\x00GFA\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00"),
-		want: NewVecDense(0, []float64{}),
+		want: &VecDense{},
 		err:  ErrZeroLength,
 		eq:   Equal,
 	},
diff --git a/mat/list_test.go b/mat/list_test.go
index 1e6e48ff..ea9c7bc2 100644
--- a/mat/list_test.go
+++ b/mat/list_test.go
@@ -315,7 +315,10 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 	case Untransposer:
 		rMatrix = retranspose(a, makeRandOf(t.Untranspose(), n, m))
 	case *Dense, *basicMatrix:
-		mat := NewDense(m, n, nil)
+		var mat = &Dense{}
+		if m != 0 && n != 0 {
+			mat = NewDense(m, n, nil)
+		}
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
 				mat.Set(i, j, rand.NormFloat64())
@@ -363,7 +366,10 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 		if m != n {
 			panic("bad size")
 		}
-		mat := NewSymDense(n, nil)
+		mat := &SymDense{}
+		if n != 0 {
+			mat = NewSymDense(n, nil)
+		}
 		for i := 0; i < m; i++ {
 			for j := i; j < n; j++ {
 				mat.SetSym(i, j, rand.NormFloat64())
@@ -386,6 +392,14 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 			triKind = (*TriDense)(t).triKind()
 		}
 
+		if n == 0 {
+			uplo := blas.Upper
+			if triKind == Lower {
+				uplo = blas.Lower
+			}
+			return returnAs(&TriDense{mat: blas64.Triangular{Uplo: uplo}}, t)
+		}
+
 		mat := NewTriDense(n, triKind, nil)
 		if triKind == Upper {
 			for i := 0; i < m; i++ {
@@ -604,7 +618,7 @@ var testMatrices = []Matrix{
 	&SymDense{},
 	NewTriDense(3, true, nil),
 	NewTriDense(3, false, nil),
-	NewVecDense(0, nil),
+	&VecDense{mat: blas64.Vector{Inc: 1}},
 	&DiagDense{},
 	&basicVector{},
 	&VecDense{mat: blas64.Vector{Inc: 10}},
@@ -618,7 +632,7 @@ var testMatrices = []Matrix{
 	TransposeTri{NewTriDense(3, true, nil)},
 	Transpose{NewTriDense(3, false, nil)},
 	TransposeTri{NewTriDense(3, false, nil)},
-	Transpose{NewVecDense(0, nil)},
+	Transpose{&VecDense{mat: blas64.Vector{Inc: 1}}},
 	Transpose{&VecDense{mat: blas64.Vector{Inc: 10}}},
 	Transpose{&DiagDense{}},
 	TransposeTri{&DiagDense{}},
diff --git a/mat/mul_test.go b/mat/mul_test.go
index c2125412..9cc9c5e4 100644
--- a/mat/mul_test.go
+++ b/mat/mul_test.go
@@ -97,7 +97,7 @@ func TestMulTypes(t *testing.T) {
 
 		// Check that it panics if it is supposed to
 		if test.Panics {
-			c := NewDense(0, 0, nil)
+			c := &Dense{}
 			fn := func() {
 				c.Mul(a, b)
 			}
@@ -125,7 +125,7 @@ func TestMulTypes(t *testing.T) {
 		ccomp := matComp{r: ar, c: bc, data: cvecCopy}
 
 		// Do normal multiply with empty dense
-		d := NewDense(0, 0, nil)
+		d := &Dense{}
 
 		testMul(t, a, b, d, acomp, bcomp, ccomp, false, "zero receiver")
 
diff --git a/mat/product_test.go b/mat/product_test.go
index 7033381b..02aae9a0 100644
--- a/mat/product_test.go
+++ b/mat/product_test.go
@@ -126,7 +126,10 @@ func TestProduct(t *testing.T) {
 			}
 		}
 
-		got := NewDense(test.product.r, test.product.c, nil)
+		got := &Dense{}
+		if test.product.r != 0 && test.product.c != 0 {
+			got = NewDense(test.product.r, test.product.c, nil)
+		}
 		panicked, message := panics(func() {
 			got.Product(factors...)
 		})
diff --git a/mat/symband.go b/mat/symband.go
index 8763a1bc..40af5bb3 100644
--- a/mat/symband.go
+++ b/mat/symband.go
@@ -67,7 +67,10 @@ type RawSymBander interface {
 // which is passed to NewBandDense as []float64{1, 2, 3, 4, ...} with k=2.
 // Only the values in the band portion of the matrix are used.
 func NewSymBandDense(n, k int, data []float64) *SymBandDense {
-	if n < 0 || k < 0 {
+	if n <= 0 || k < 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
 		panic("mat: negative dimension")
 	}
 	if k+1 > n {
diff --git a/mat/symmetric.go b/mat/symmetric.go
index d3ac2617..2e2fc1c7 100644
--- a/mat/symmetric.go
+++ b/mat/symmetric.go
@@ -55,12 +55,16 @@ type MutableSymmetric interface {
 // a new slice is allocated for the backing slice. If len(data) == n*n, data is
 // used as the backing slice, and changes to the elements of the returned SymDense
 // will be reflected in data. If neither of these is true, NewSymDense will panic.
+// NewSymDense will panic if n is zero.
 //
 // The data must be arranged in row-major order, i.e. the (i*c + j)-th
 // element in the data slice is the {i, j}-th element in the matrix.
 // Only the values in the upper triangular portion of the matrix are used.
 func NewSymDense(n int, data []float64) *SymDense {
-	if n < 0 {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
 		panic("mat: negative dimension")
 	}
 	if data != nil && n*n != len(data) {
diff --git a/mat/triangular.go b/mat/triangular.go
index c4446b55..751a58e8 100644
--- a/mat/triangular.go
+++ b/mat/triangular.go
@@ -111,12 +111,16 @@ func (t TransposeTri) UntransposeTri() Triangular {
 // a new slice is allocated for the backing slice. If len(data) == n*n, data is
 // used as the backing slice, and changes to the elements of the returned TriDense
 // will be reflected in data. If neither of these is true, NewTriDense will panic.
+// NewTriDense will panic if n is zero.
 //
 // The data must be arranged in row-major order, i.e. the (i*c + j)-th
 // element in the data slice is the {i, j}-th element in the matrix.
 // Only the values in the triangular portion corresponding to kind are used.
 func NewTriDense(n int, kind TriKind, data []float64) *TriDense {
-	if n < 0 {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
 		panic("mat: negative dimension")
 	}
 	if data != nil && len(data) != n*n {
diff --git a/mat/vector.go b/mat/vector.go
index ca29e503..d2409131 100644
--- a/mat/vector.go
+++ b/mat/vector.go
@@ -86,8 +86,12 @@ type VecDense struct {
 // a new slice is allocated for the backing slice. If len(data) == n, data is
 // used as the backing slice, and changes to the elements of the returned VecDense
 // will be reflected in data. If neither of these is true, NewVecDense will panic.
+// NewVecDense will panic if n is zero.
 func NewVecDense(n int, data []float64) *VecDense {
-	if n < 0 {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
 		panic("mat: negative dimension")
 	}
 	if len(data) != n && data != nil {