mat: Expose ReuseAs method to help resizing (#1082)

* mat: Expose ReuseAs method to help resizing

Fixes #963

mat/band.go

@@ -10,9 +10,11 @@ import (
 
 var (
 	bandDense *BandDense
-	_         Matrix    = bandDense
-	_         Banded    = bandDense
-	_         RawBander = bandDense
+	_         Matrix      = bandDense
+	_         allMatrix   = bandDense
+	_         denseMatrix = bandDense
+	_         Banded      = bandDense
+	_         RawBander   = bandDense
 
 	_ NonZeroDoer    = bandDense
 	_ RowNonZeroDoer = bandDense
@@ -195,6 +197,26 @@ func (b *BandDense) SetRawBand(mat blas64.Band) {
 	b.mat = mat
 }
 
+// IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
+// receiver for size-restricted operations. Dense matrices can be zeroed using Reset.
+func (b *BandDense) IsZero() bool {
+	return b.mat.Stride == 0
+}
+
+// Reset zeros the dimensions of the matrix so that it can be reused as the
+// receiver of a dimensionally restricted operation.
+//
+// Reset should not be used when the matrix shares backing data.
+// See the Reseter interface for more information.
+func (b *BandDense) Reset() {
+	b.mat.Rows = 0
+	b.mat.Cols = 0
+	b.mat.KL = 0
+	b.mat.KU = 0
+	b.mat.Stride = 0
+	b.mat.Data = b.mat.Data[:0:0]
+}
+
 // DiagView returns the diagonal as a matrix backed by the original data.
 func (b *BandDense) DiagView() Diagonal {
 	n := min(b.mat.Rows, b.mat.Cols)

mat/cdense.go

@@ -6,7 +6,14 @@ package mat
 
 import "gonum.org/v1/gonum/blas/cblas128"
 
-// Dense is a dense matrix representation with complex data.
+var (
+	cDense *CDense
+
+	_ CMatrix   = cDense
+	_ allMatrix = cDense
+)
+
+// CDense is a dense matrix representation with complex data.
 type CDense struct {
 	mat cblas128.General
 
@@ -62,7 +69,7 @@ func NewCDense(r, c int, data []complex128) *CDense {
 // or checks that a non-empty matrix is r×c.
 //
 // reuseAs must be kept in sync with reuseAsZeroed.
-func (m *CDense) reuseAs(r, c int) {
+func (m *CDense) reuseAsNonZeroed(r, c int) {
 	if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
 		// Panic as a string, not a mat.Error.
 		panic("mat: caps not correctly set")
@@ -115,6 +122,7 @@ func (m *CDense) reuseAsZeroed(r, c int) {
 // Reset zeros the dimensions of the matrix so that it can be reused as the
 // receiver of a dimensionally restricted operation.
 //
+// Reset should not be used when the matrix shares backing data.
 // See the Reseter interface for more information.
 func (m *CDense) Reset() {
 	// Row, Cols and Stride must be zeroed in unison.

mat/cholesky.go

@@ -220,7 +220,7 @@ func (c *Cholesky) SolveTo(dst *Dense, b Matrix) error {
 		panic(ErrShape)
 	}
 
-	dst.reuseAs(bm, bn)
+	dst.reuseAsNonZeroed(bm, bn)
 	if b != dst {
 		dst.Copy(b)
 	}
@@ -267,14 +267,14 @@ func (c *Cholesky) SolveVecTo(dst *VecDense, b Vector) error {
 	}
 	switch rv := b.(type) {
 	default:
-		dst.reuseAs(n)
+		dst.reuseAsNonZeroed(n)
 		return c.SolveTo(dst.asDense(), b)
 	case RawVectorer:
 		bmat := rv.RawVector()
 		if dst != b {
 			dst.checkOverlap(bmat)
 		}
-		dst.reuseAs(n)
+		dst.reuseAsNonZeroed(n)
 		if dst != b {
 			dst.CopyVec(b)
 		}
@@ -305,7 +305,7 @@ func (c *Cholesky) UTo(dst *TriDense) *TriDense {
 	if dst == nil {
 		dst = NewTriDense(n, Upper, make([]float64, n*n))
 	} else {
-		dst.reuseAs(n, Upper)
+		dst.reuseAsNonZeroed(n, Upper)
 	}
 	dst.Copy(c.chol)
 	return dst
@@ -323,7 +323,7 @@ func (c *Cholesky) LTo(dst *TriDense) *TriDense {
 	if dst == nil {
 		dst = NewTriDense(n, Lower, make([]float64, n*n))
 	} else {
-		dst.reuseAs(n, Lower)
+		dst.reuseAsNonZeroed(n, Lower)
 	}
 	dst.Copy(c.chol.TTri())
 	return dst
@@ -340,7 +340,7 @@ func (c *Cholesky) ToSym(dst *SymDense) *SymDense {
 	if dst == nil {
 		dst = NewSymDense(n, nil)
 	} else {
-		dst.reuseAs(n)
+		dst.reuseAsNonZeroed(n)
 	}
 	// Create a TriDense representing the Cholesky factor U with dst's
 	// backing slice.
@@ -378,7 +378,7 @@ func (c *Cholesky) InverseTo(s *SymDense) error {
 	if !c.valid() {
 		panic(badCholesky)
 	}
-	s.reuseAs(c.chol.mat.N)
+	s.reuseAsNonZeroed(c.chol.mat.N)
 	// Create a TriDense representing the Cholesky factor U with the backing
 	// slice from s.
 	// Operations on u are reflected in s.

mat/dense.go

@@ -12,8 +12,10 @@ import (
 var (
 	dense *Dense
 
-	_ Matrix  = dense
-	_ Mutable = dense
+	_ Matrix      = dense
+	_ allMatrix   = dense
+	_ denseMatrix = dense
+	_ Mutable     = dense
 
 	_ ClonerFrom   = dense
 	_ RowViewer    = dense
@@ -47,7 +49,7 @@ func NewDense(r, c int, data []float64) *Dense {
 		if r == 0 || c == 0 {
 			panic(ErrZeroLength)
 		}
-		panic("mat: negative dimension")
+		panic(ErrNegativeDimension)
 	}
 	if data != nil && r*c != len(data) {
 		panic(ErrShape)
@@ -67,11 +69,32 @@ func NewDense(r, c int, data []float64) *Dense {
 	}
 }
 
-// reuseAs resizes an empty matrix to a r×c matrix,
-// or checks that a non-empty matrix is r×c.
+// ReuseAs changes the receiver if it IsZero() to be of size r×c.
 //
-// reuseAs must be kept in sync with reuseAsZeroed.
-func (m *Dense) reuseAs(r, c int) {
+// ReuseAs re-uses the backing data slice if it has sufficient capacity,
+// otherwise a new slice is allocated. The data is then zeroed.
+//
+// ReuseAs panics if the receiver is not zero-sized, and panics if
+// the input sizes are less than one. To zero the receiver for re-use,
+// Reset should be used.
+func (m *Dense) ReuseAs(r, c int) {
+	if r <= 0 || c <= 0 {
+		if r == 0 || c == 0 {
+			panic(ErrZeroLength)
+		}
+		panic(ErrNegativeDimension)
+	}
+	if !m.IsZero() {
+		panic(ErrReuseNonZero)
+	}
+	m.reuseAsZeroed(r, c)
+}
+
+// reuseAsNonZeroed resizes an empty matrix to a r×c matrix,
+// or checks that a non-empty matrix is r×c. It does not zero
+// the data in the receiver.
+func (m *Dense) reuseAsNonZeroed(r, c int) {
+	// reuseAs must be kept in sync with reuseAsZeroed.
 	if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
 		// Panic as a string, not a mat.Error.
 		panic("mat: caps not correctly set")
@@ -98,9 +121,8 @@ func (m *Dense) reuseAs(r, c int) {
 // reuseAsZeroed resizes an empty matrix to a r×c matrix,
 // or checks that a non-empty matrix is r×c. It zeroes
 // all the elements of the matrix.
-//
-// reuseAsZeroed must be kept in sync with reuseAs.
 func (m *Dense) reuseAsZeroed(r, c int) {
+	// reuseAsZeroed must be kept in sync with reuseAsNonZeroed.
 	if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
 		// Panic as a string, not a mat.Error.
 		panic("mat: caps not correctly set")
@@ -152,6 +174,7 @@ func (m *Dense) isolatedWorkspace(a Matrix) (w *Dense, restore func()) {
 // Reset zeros the dimensions of the matrix so that it can be reused as the
 // receiver of a dimensionally restricted operation.
 //
+// Reset should not be used when the matrix shares backing data.
 // See the Reseter interface for more information.
 func (m *Dense) Reset() {
 	// Row, Cols and Stride must be zeroed in unison.
@@ -518,7 +541,7 @@ func (m *Dense) Stack(a, b Matrix) {
 		panic(ErrShape)
 	}
 
-	m.reuseAs(ar+br, ac)
+	m.reuseAsNonZeroed(ar+br, ac)
 
 	m.Copy(a)
 	w := m.Slice(ar, ar+br, 0, bc).(*Dense)
@@ -536,7 +559,7 @@ func (m *Dense) Augment(a, b Matrix) {
 		panic(ErrShape)
 	}
 
-	m.reuseAs(ar, ac+bc)
+	m.reuseAsNonZeroed(ar, ac+bc)
 
 	m.Copy(a)
 	w := m.Slice(0, br, ac, ac+bc).(*Dense)

mat/dense_arithmetic.go

@@ -23,7 +23,7 @@ func (m *Dense) Add(a, b Matrix) {
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
-	m.reuseAs(ar, ac)
+	m.reuseAsNonZeroed(ar, ac)
 
 	if arm, ok := a.(*Dense); ok {
 		if brm, ok := b.(*Dense); ok {
@@ -72,7 +72,7 @@ func (m *Dense) Sub(a, b Matrix) {
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
-	m.reuseAs(ar, ac)
+	m.reuseAsNonZeroed(ar, ac)
 
 	if arm, ok := a.(*Dense); ok {
 		if brm, ok := b.(*Dense); ok {
@@ -122,7 +122,7 @@ func (m *Dense) MulElem(a, b Matrix) {
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
-	m.reuseAs(ar, ac)
+	m.reuseAsNonZeroed(ar, ac)
 
 	if arm, ok := a.(*Dense); ok {
 		if brm, ok := b.(*Dense); ok {
@@ -172,7 +172,7 @@ func (m *Dense) DivElem(a, b Matrix) {
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
-	m.reuseAs(ar, ac)
+	m.reuseAsNonZeroed(ar, ac)
 
 	if arm, ok := a.(*Dense); ok {
 		if brm, ok := b.(*Dense); ok {
@@ -220,7 +220,7 @@ func (m *Dense) Inverse(a Matrix) error {
 	if r != c {
 		panic(ErrSquare)
 	}
-	m.reuseAs(a.Dims())
+	m.reuseAsNonZeroed(a.Dims())
 	aU, aTrans := untransposeExtract(a)
 	switch rm := aU.(type) {
 	case *Dense:
@@ -278,7 +278,7 @@ func (m *Dense) Mul(a, b Matrix) {
 
 	aU, aTrans := untransposeExtract(a)
 	bU, bTrans := untransposeExtract(b)
-	m.reuseAs(ar, bc)
+	m.reuseAsNonZeroed(ar, bc)
 	var restore func()
 	if m == aU {
 		m, restore = m.isolatedWorkspace(aU)
@@ -471,7 +471,7 @@ func (m *Dense) Exp(a Matrix) {
 		panic(ErrShape)
 	}
 
-	m.reuseAs(r, r)
+	m.reuseAsNonZeroed(r, r)
 	if r == 1 {
 		m.mat.Data[0] = math.Exp(a.At(0, 0))
 		return
@@ -633,7 +633,7 @@ func (m *Dense) Pow(a Matrix, n int) {
 		panic(ErrShape)
 	}
 
-	m.reuseAs(r, c)
+	m.reuseAsNonZeroed(r, c)
 
 	// Take possible fast paths.
 	switch n {
@@ -679,7 +679,7 @@ func (m *Dense) Pow(a Matrix, n int) {
 func (m *Dense) Scale(f float64, a Matrix) {
 	ar, ac := a.Dims()
 
-	m.reuseAs(ar, ac)
+	m.reuseAsNonZeroed(ar, ac)
 
 	aU, aTrans := untransposeExtract(a)
 	if rm, ok := aU.(*Dense); ok {
@@ -719,7 +719,7 @@ func (m *Dense) Scale(f float64, a Matrix) {
 func (m *Dense) Apply(fn func(i, j int, v float64) float64, a Matrix) {
 	ar, ac := a.Dims()
 
-	m.reuseAs(ar, ac)
+	m.reuseAsNonZeroed(ar, ac)
 
 	aU, aTrans := untransposeExtract(a)
 	if rm, ok := aU.(*Dense); ok {
@@ -794,14 +794,14 @@ func (m *Dense) RankOne(a Matrix, alpha float64, x, y Vector) {
 
 	if fast {
 		if m != a {
-			m.reuseAs(ar, ac)
+			m.reuseAsNonZeroed(ar, ac)
 			m.Copy(a)
 		}
 		blas64.Ger(alpha, xmat, ymat, m.mat)
 		return
 	}
 
-	m.reuseAs(ar, ac)
+	m.reuseAsNonZeroed(ar, ac)
 	for i := 0; i < ar; i++ {
 		for j := 0; j < ac; j++ {
 			m.set(i, j, a.At(i, j)+alpha*x.AtVec(i)*y.AtVec(j))

mat/diagonal.go

@@ -12,6 +12,8 @@ import (
 var (
 	diagDense *DiagDense
 	_         Matrix          = diagDense
+	_         allMatrix       = diagDense
+	_         denseMatrix     = diagDense
 	_         Diagonal        = diagDense
 	_         MutableDiagonal = diagDense
 	_         Triangular      = diagDense
@@ -156,9 +158,10 @@ func (d *DiagDense) TriBand() (n, k int, kind TriKind) {
 	return d.mat.N, 0, Upper
 }
 
-// Reset zeros the length of the matrix so that it can be reused as the
+// Reset zeros the dimensions of the matrix so that it can be reused as the
 // receiver of a dimensionally restricted operation.
 //
+// Reset should not be used when the matrix shares backing data.
 // See the Reseter interface for more information.
 func (d *DiagDense) Reset() {
 	// No change of Inc or n to 0 may be
@@ -184,7 +187,7 @@ func (d *DiagDense) DiagView() Diagonal {
 // be min(r, c) long or zero. Otherwise DiagFrom will panic.
 func (d *DiagDense) DiagFrom(m Matrix) {
 	n := min(m.Dims())
-	d.reuseAs(n)
+	d.reuseAsNonZeroed(n)
 
 	var vec blas64.Vector
 	switch r := m.(type) {
@@ -283,9 +286,9 @@ func (d *DiagDense) RawSymBand() blas64.SymmetricBand {
 	}
 }
 
-// reuseAs resizes an empty diagonal to a r×r diagonal,
+// reuseAsNonZeroed resizes an empty diagonal to a r×r diagonal,
 // or checks that a non-empty matrix is r×r.
-func (d *DiagDense) reuseAs(r int) {
+func (d *DiagDense) reuseAsNonZeroed(r int) {
 	if r == 0 {
 		panic(ErrZeroLength)
 	}

mat/doc.go

@@ -51,6 +51,8 @@
 //  - var a mat.Dense
 //  - a := NewDense(0, 0, make([]float64, 0, 100))
 //  - a.Reset()
+// Reset should not be used when multiple different matrices share the same backing
+// data slice. This can cause unexpected data modifications after being resized.
 // A zero-value matrix can not be sliced even if it does have an adequately sized
 // backing data slice, but can be expanded using its Grow method if it exists.
 //

mat/eigen.go

@@ -109,7 +109,7 @@ func (e *EigenSym) VectorsTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(r, c)
+		dst.reuseAsNonZeroed(r, c)
 	}
 	dst.Copy(e.vectors)
 	return dst
@@ -321,7 +321,7 @@ func (e *Eigen) VectorsTo(dst *CDense) *CDense {
 	if dst == nil {
 		dst = NewCDense(e.n, e.n, nil)
 	} else {
-		dst.reuseAs(e.n, e.n)
+		dst.reuseAsNonZeroed(e.n, e.n)
 	}
 	dst.Copy(e.rVectors)
 	return dst
@@ -343,7 +343,7 @@ func (e *Eigen) LeftVectorsTo(dst *CDense) *CDense {
 	if dst == nil {
 		dst = NewCDense(e.n, e.n, nil)
 	} else {
-		dst.reuseAs(e.n, e.n)
+		dst.reuseAsNonZeroed(e.n, e.n)
 	}
 	dst.Copy(e.lVectors)
 	return dst

mat/errors.go

@@ -114,7 +114,9 @@ type Error struct{ string }
 func (err Error) Error() string { return err.string }
 
 var (
+	ErrNegativeDimension   = Error{"matrix: negative dimension"}
 	ErrIndexOutOfRange     = Error{"matrix: index out of range"}
+	ErrReuseNonZero        = Error{"mat: reuse of non-zero matrix"}
 	ErrRowAccess           = Error{"matrix: row index out of range"}
 	ErrColAccess           = Error{"matrix: column index out of range"}
 	ErrVectorAccess        = Error{"matrix: vector index out of range"}

mat/gsvd.go

@@ -344,7 +344,7 @@ func (gsvd *GSVD) UTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(r, c)
+		dst.reuseAsNonZeroed(r, c)
 	}
 
 	tmp := &Dense{
@@ -373,7 +373,7 @@ func (gsvd *GSVD) VTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(r, c)
+		dst.reuseAsNonZeroed(r, c)
 	}
 
 	tmp := &Dense{
@@ -402,7 +402,7 @@ func (gsvd *GSVD) QTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(r, c)
+		dst.reuseAsNonZeroed(r, c)
 	}
 
 	tmp := &Dense{

mat/hogsvd.go

@@ -173,7 +173,7 @@ func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense {
 		r, c := gsvd.b[n].Dims()
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(gsvd.b[n].Dims())
+		dst.reuseAsNonZeroed(gsvd.b[n].Dims())
 	}
 	dst.Copy(&gsvd.b[n])
 	var v VecDense
@@ -226,7 +226,7 @@ func (gsvd *HOGSVD) VTo(dst *Dense) *Dense {
 		r, c := gsvd.v.Dims()
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(gsvd.v.Dims())
+		dst.reuseAsNonZeroed(gsvd.v.Dims())
 	}
 	dst.Copy(gsvd.v)
 	return dst

mat/io.go

@@ -175,7 +175,7 @@ func (m *Dense) UnmarshalBinary(data []byte) error {
 	}
 
 	p := headerSize
-	m.reuseAs(int(rows), int(cols))
+	m.reuseAsNonZeroed(int(rows), int(cols))
 	for i := range m.mat.Data {
 		m.mat.Data[i] = math.Float64frombits(binary.LittleEndian.Uint64(data[p : p+sizeFloat64]))
 		p += sizeFloat64
@@ -226,7 +226,7 @@ func (m *Dense) UnmarshalBinaryFrom(r io.Reader) (int, error) {
 		return n, errTooBig
 	}
 
-	m.reuseAs(int(rows), int(cols))
+	m.reuseAsNonZeroed(int(rows), int(cols))
 	var b [8]byte
 	for i := range m.mat.Data {
 		nn, err := readFull(r, b[:])
@@ -362,7 +362,7 @@ func (v *VecDense) UnmarshalBinary(data []byte) error {
 	}
 
 	p := headerSize
-	v.reuseAs(int(n))
+	v.reuseAsNonZeroed(int(n))
 	for i := range v.mat.Data {
 		v.mat.Data[i] = math.Float64frombits(binary.LittleEndian.Uint64(data[p : p+sizeFloat64]))
 		p += sizeFloat64
@@ -407,7 +407,7 @@ func (v *VecDense) UnmarshalBinaryFrom(r io.Reader) (int, error) {
 		return n, errTooBig
 	}
 
-	v.reuseAs(int(l))
+	v.reuseAsNonZeroed(int(l))
 	var b [8]byte
 	for i := range v.mat.Data {
 		nn, err := readFull(r, b[:])

mat/lq.go

@@ -98,7 +98,7 @@ func (lq *LQ) LTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(r, c)
+		dst.reuseAsNonZeroed(r, c)
 	}
 
 	// Disguise the LQ as a lower triangular.
@@ -183,12 +183,12 @@ func (lq *LQ) SolveTo(dst *Dense, trans bool, b Matrix) error {
 		if c != br {
 			panic(ErrShape)
 		}
-		dst.reuseAs(r, bc)
+		dst.reuseAsNonZeroed(r, bc)
 	} else {
 		if r != br {
 			panic(ErrShape)
 		}
-		dst.reuseAs(c, bc)
+		dst.reuseAsNonZeroed(c, bc)
 	}
 	// Do not need to worry about overlap between x and b because w has its own
 	// independent storage.
@@ -254,9 +254,9 @@ func (lq *LQ) SolveVecTo(dst *VecDense, trans bool, b Vector) error {
 		bm = b.asDense()
 	}
 	if trans {
-		dst.reuseAs(r)
+		dst.reuseAsNonZeroed(r)
 	} else {
-		dst.reuseAs(c)
+		dst.reuseAsNonZeroed(c)
 	}
 	return lq.SolveTo(dst.asDense(), trans, bm)
 }

mat/lu.go

@@ -72,7 +72,7 @@ func (lu *LU) factorize(a Matrix, norm lapack.MatrixNorm) {
 		lu.lu = NewDense(r, r, nil)
 	} else {
 		lu.lu.Reset()
-		lu.lu.reuseAs(r, r)
+		lu.lu.reuseAsNonZeroed(r, r)
 	}
 	lu.lu.Copy(a)
 	if cap(lu.pivot) < r {
@@ -205,7 +205,7 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y Vector) {
 			if lu.lu == nil {
 				lu.lu = NewDense(n, n, nil)
 			} else {
-				lu.lu.reuseAs(n, n)
+				lu.lu.reuseAsNonZeroed(n, n)
 			}
 		} else if len(lu.pivot) != n {
 			panic(ErrShape)
@@ -261,7 +261,7 @@ func (lu *LU) LTo(dst *TriDense) *TriDense {
 	if dst == nil {
 		dst = NewTriDense(n, Lower, nil)
 	} else {
-		dst.reuseAs(n, Lower)
+		dst.reuseAsNonZeroed(n, Lower)
 	}
 	// Extract the lower triangular elements.
 	for i := 0; i < n; i++ {
@@ -288,7 +288,7 @@ func (lu *LU) UTo(dst *TriDense) *TriDense {
 	if dst == nil {
 		dst = NewTriDense(n, Upper, nil)
 	} else {
-		dst.reuseAs(n, Upper)
+		dst.reuseAsNonZeroed(n, Upper)
 	}
 	// Extract the upper triangular elements.
 	for i := 0; i < n; i++ {
@@ -304,7 +304,7 @@ func (lu *LU) UTo(dst *TriDense) *TriDense {
 // and all other elements equal to zero. swaps[i] specifies the row with which
 // i will be swapped, which is equivalent to the non-zero column of row i.
 func (m *Dense) Permutation(r int, swaps []int) {
-	m.reuseAs(r, r)
+	m.reuseAsNonZeroed(r, r)
 	for i := 0; i < r; i++ {
 		zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+r])
 		v := swaps[i]
@@ -341,7 +341,7 @@ func (lu *LU) SolveTo(dst *Dense, trans bool, b Matrix) error {
 		return Condition(math.Inf(1))
 	}
 
-	dst.reuseAs(n, bc)
+	dst.reuseAsNonZeroed(n, bc)
 	bU, _ := untranspose(b)
 	var restore func()
 	if dst == bU {
@@ -384,7 +384,7 @@ func (lu *LU) SolveVecTo(dst *VecDense, trans bool, b Vector) error {
 	}
 	switch rv := b.(type) {
 	default:
-		dst.reuseAs(n)
+		dst.reuseAsNonZeroed(n)
 		return lu.SolveTo(dst.asDense(), trans, b)
 	case RawVectorer:
 		if dst != b {
@@ -396,7 +396,7 @@ func (lu *LU) SolveVecTo(dst *VecDense, trans bool, b Vector) error {
 			return Condition(math.Inf(1))
 		}
 
-		dst.reuseAs(n)
+		dst.reuseAsNonZeroed(n)
 		var restore func()
 		if dst == b {
 			dst, restore = dst.isolatedWorkspace(b)

mat/matrix.go

@@ -30,6 +30,22 @@ type Matrix interface {
 	T() Matrix
 }
 
+// allMatrix represents the extra set of methods that all mat Matrix types
+// should satisfy. This is used to enforce compile-time consistency between the
+// Dense types, especially helpful when adding new features.
+type allMatrix interface {
+	Reseter
+	IsZero() bool
+}
+
+// denseMatrix represents the extra set of methods that all Dense Matrix types
+// should satisfy. This is used to enforce compile-time consistency between the
+// Dense types, especially helpful when adding new features.
+type denseMatrix interface {
+	DiagView() Diagonal
+	Tracer
+}
+
 var (
 	_ Matrix       = Transpose{}
 	_ Untransposer = Transpose{}
@@ -140,8 +156,9 @@ type ClonerFrom interface {
 // restricted operation. This is commonly used when the matrix is being used as a workspace
 // or temporary matrix.
 //
-// If the matrix is a view, using the reset matrix may result in data corruption in elements
-// outside the view.
+// If the matrix is a view, using Reset may result in data corruption in elements outside
+// the view. Similarly, if the matrix shares backing data with another variable, using
+// Reset may lead to unexpected changes in data values.
 type Reseter interface {
 	Reset()
 }

mat/product.go

@@ -31,7 +31,7 @@ func (m *Dense) Product(factors ...Matrix) {
 		}
 		return
 	case 1:
-		m.reuseAs(factors[0].Dims())
+		m.reuseAsNonZeroed(factors[0].Dims())
 		m.Copy(factors[0])
 		return
 	case 2:
@@ -43,7 +43,7 @@ func (m *Dense) Product(factors ...Matrix) {
 	p := newMultiplier(m, factors)
 	p.optimize()
 	result := p.multiply()
-	m.reuseAs(result.Dims())
+	m.reuseAsNonZeroed(result.Dims())
 	m.Copy(result)
 	putWorkspace(result)
 }

mat/qr.go

@@ -98,7 +98,7 @@ func (qr *QR) RTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(r, c)
+		dst.reuseAsNonZeroed(r, c)
 	}
 
 	// Disguise the QR as an upper triangular
@@ -178,12 +178,12 @@ func (qr *QR) SolveTo(dst *Dense, trans bool, b Matrix) error {
 		if c != br {
 			panic(ErrShape)
 		}
-		dst.reuseAs(r, bc)
+		dst.reuseAsNonZeroed(r, bc)
 	} else {
 		if r != br {
 			panic(ErrShape)
 		}
-		dst.reuseAs(c, bc)
+		dst.reuseAsNonZeroed(c, bc)
 	}
 	// Do not need to worry about overlap between m and b because x has its own
 	// independent storage.
@@ -250,9 +250,9 @@ func (qr *QR) SolveVecTo(dst *VecDense, trans bool, b Vector) error {
 		bm = b.asDense()
 	}
 	if trans {
-		dst.reuseAs(r)
+		dst.reuseAsNonZeroed(r)
 	} else {
-		dst.reuseAs(c)
+		dst.reuseAsNonZeroed(c)
 	}
 	return qr.SolveTo(dst.asDense(), trans, bm)
 

mat/solve.go

@@ -24,7 +24,7 @@ func (m *Dense) Solve(a, b Matrix) error {
 	if ar != br {
 		panic(ErrShape)
 	}
-	m.reuseAs(ac, bc)
+	m.reuseAsNonZeroed(ac, bc)
 
 	// TODO(btracey): Add special cases for SymDense, etc.
 	aU, aTrans := untranspose(a)
@@ -121,7 +121,7 @@ func (v *VecDense) SolveVec(a Matrix, b Vector) error {
 		if v != b {
 			v.checkOverlap(bmat)
 		}
-		v.reuseAs(c)
+		v.reuseAsNonZeroed(c)
 		m := v.asDense()
 		// We conditionally create bm as m when b and v are identical
 		// to prevent the overlap detection code from identifying m
@@ -134,7 +134,7 @@ func (v *VecDense) SolveVec(a Matrix, b Vector) error {
 		return m.Solve(a, bm)
 	}
 
-	v.reuseAs(c)
+	v.reuseAsNonZeroed(c)
 	m := v.asDense()
 	return m.Solve(a, b)
 }

mat/svd.go

@@ -199,7 +199,7 @@ func (svd *SVD) UTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(r, c, nil)
 	} else {
-		dst.reuseAs(r, c)
+		dst.reuseAsNonZeroed(r, c)
 	}
 
 	tmp := &Dense{
@@ -233,7 +233,7 @@ func (svd *SVD) VTo(dst *Dense) *Dense {
 	if dst == nil {
 		dst = NewDense(c, r, nil)
 	} else {
-		dst.reuseAs(c, r)
+		dst.reuseAsNonZeroed(c, r)
 	}
 
 	tmp := &Dense{

mat/symband.go

@@ -12,6 +12,8 @@ import (
 var (
 	symBandDense *SymBandDense
 	_            Matrix           = symBandDense
+	_            allMatrix        = symBandDense
+	_            denseMatrix      = symBandDense
 	_            Symmetric        = symBandDense
 	_            Banded           = symBandDense
 	_            SymBanded        = symBandDense
@@ -157,6 +159,25 @@ func (s *SymBandDense) SetRawSymBand(mat blas64.SymmetricBand) {
 	s.mat = mat
 }
 
+// IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
+// receiver for size-restricted operations. Dense matrices can be zeroed using Reset.
+func (s *SymBandDense) IsZero() bool {
+	return s.mat.Stride == 0
+}
+
+// Reset zeros the dimensions of the matrix so that it can be reused as the
+// receiver of a dimensionally restricted operation.
+//
+// Reset should not be used when the matrix shares backing data.
+// See the Reseter interface for more information.
+func (s *SymBandDense) Reset() {
+	s.mat.N = 0
+	s.mat.K = 0
+	s.mat.Stride = 0
+	s.mat.Uplo = 0
+	s.mat.Data = s.mat.Data[:0:0]
+}
+
 // Zero sets all of the matrix elements to zero.
 func (s *SymBandDense) Zero() {
 	for i := 0; i < s.mat.N; i++ {

mat/symmetric.go

@@ -15,6 +15,8 @@ var (
 	symDense *SymDense
 
 	_ Matrix           = symDense
+	_ allMatrix        = symDense
+	_ denseMatrix      = symDense
 	_ Symmetric        = symDense
 	_ RawSymmetricer   = symDense
 	_ MutableSymmetric = symDense
@@ -127,6 +129,7 @@ func (s *SymDense) SetRawSymmetric(mat blas64.Symmetric) {
 // Reset zeros the dimensions of the matrix so that it can be reused as the
 // receiver of a dimensionally restricted operation.
 //
+// Reset should not be used when the matrix shares backing data.
 // See the Reseter interface for more information.
 func (s *SymDense) Reset() {
 	// N and Stride must be zeroed in unison.
@@ -134,6 +137,27 @@ func (s *SymDense) Reset() {
 	s.mat.Data = s.mat.Data[:0]
 }
 
+// ReuseAsSym changes the receiver if it IsZero() to be of size n×n.
+//
+// ReuseAsSym re-uses the backing data slice if it has sufficient capacity,
+// otherwise a new slice is allocated. The data is then zeroed.
+//
+// ReuseAsSym panics if the receiver is not zero-sized, and panics if
+// the input size is less than one. To zero the receiver for re-use,
+// Reset should be used.
+func (s *SymDense) ReuseAsSym(n int) {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
+		panic(ErrNegativeDimension)
+	}
+	if !s.IsZero() {
+		panic(ErrReuseNonZero)
+	}
+	s.reuseAsZeroed(n)
+}
+
 // Zero sets all of the matrix elements to zero.
 func (s *SymDense) Zero() {
 	for i := 0; i < s.mat.N; i++ {
@@ -149,9 +173,10 @@ func (s *SymDense) IsZero() bool {
 	return s.mat.N == 0
 }
 
-// reuseAs resizes an empty matrix to a n×n matrix,
+// reuseAsNonZeroed resizes an empty matrix to a n×n matrix,
 // or checks that a non-empty matrix is n×n.
-func (s *SymDense) reuseAs(n int) {
+func (s *SymDense) reuseAsNonZeroed(n int) {
+	// reuseAsNonZeroed must be kept in sync with reuseAsZeroed.
 	if n == 0 {
 		panic(ErrZeroLength)
 	}
@@ -176,6 +201,36 @@ func (s *SymDense) reuseAs(n int) {
 	}
 }
 
+// reuseAsNonZeroed resizes an empty matrix to a n×n matrix,
+// or checks that a non-empty matrix is n×n. It then zeros the
+// elements of the matrix.
+func (s *SymDense) reuseAsZeroed(n int) {
+	// reuseAsZeroed must be kept in sync with reuseAsNonZeroed.
+	if n == 0 {
+		panic(ErrZeroLength)
+	}
+	if s.mat.N > s.cap {
+		panic(badSymCap)
+	}
+	if s.IsZero() {
+		s.mat = blas64.Symmetric{
+			N:      n,
+			Stride: n,
+			Data:   useZeroed(s.mat.Data, n*n),
+			Uplo:   blas.Upper,
+		}
+		s.cap = n
+		return
+	}
+	if s.mat.Uplo != blas.Upper {
+		panic(badSymTriangle)
+	}
+	if s.mat.N != n {
+		panic(ErrShape)
+	}
+	s.Zero()
+}
+
 func (s *SymDense) isolatedWorkspace(a Symmetric) (w *SymDense, restore func()) {
 	n := a.Symmetric()
 	if n == 0 {
@@ -205,7 +260,7 @@ func (s *SymDense) AddSym(a, b Symmetric) {
 	if n != b.Symmetric() {
 		panic(ErrShape)
 	}
-	s.reuseAs(n)
+	s.reuseAsNonZeroed(n)
 
 	if a, ok := a.(RawSymmetricer); ok {
 		if b, ok := b.(RawSymmetricer); ok {
@@ -271,7 +326,7 @@ func (s *SymDense) SymRankOne(a Symmetric, alpha float64, x Vector) {
 	if a.Symmetric() != n {
 		panic(ErrShape)
 	}
-	s.reuseAs(n)
+	s.reuseAsNonZeroed(n)
 
 	if s != a {
 		if rs, ok := a.(RawSymmetricer); ok {
@@ -317,7 +372,7 @@ func (s *SymDense) SymRankK(a Symmetric, alpha float64, x Matrix) {
 		if rs, ok := a.(RawSymmetricer); ok {
 			s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
 		}
-		s.reuseAs(n)
+		s.reuseAsNonZeroed(n)
 		s.CopySym(a)
 	}
 	t := blas.NoTrans
@@ -415,13 +470,13 @@ func (s *SymDense) RankTwo(a Symmetric, alpha float64, x, y Vector) {
 		if rs, ok := a.(RawSymmetricer); ok {
 			s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
 		}
-		s.reuseAs(n)
+		s.reuseAsNonZeroed(n)
 		s.CopySym(a)
 	}
 
 	if fast {
 		if s != a {
-			s.reuseAs(n)
+			s.reuseAsNonZeroed(n)
 			s.CopySym(a)
 		}
 		blas64.Syr2(alpha, xmat, ymat, s.mat)
@@ -429,7 +484,7 @@ func (s *SymDense) RankTwo(a Symmetric, alpha float64, x, y Vector) {
 	}
 
 	for i := 0; i < n; i++ {
-		s.reuseAs(n)
+		s.reuseAsNonZeroed(n)
 		for j := i; j < n; j++ {
 			s.set(i, j, a.At(i, j)+alpha*(x.AtVec(i)*y.AtVec(j)+y.AtVec(i)*x.AtVec(j)))
 		}
@@ -439,7 +494,7 @@ func (s *SymDense) RankTwo(a Symmetric, alpha float64, x, y Vector) {
 // ScaleSym multiplies the elements of a by f, placing the result in the receiver.
 func (s *SymDense) ScaleSym(f float64, a Symmetric) {
 	n := a.Symmetric()
-	s.reuseAs(n)
+	s.reuseAsNonZeroed(n)
 	if a, ok := a.(RawSymmetricer); ok {
 		amat := a.RawSymmetric()
 		if s != a {
@@ -467,7 +522,7 @@ func (s *SymDense) ScaleSym(f float64, a Symmetric) {
 func (s *SymDense) SubsetSym(a Symmetric, set []int) {
 	n := len(set)
 	na := a.Symmetric()
-	s.reuseAs(n)
+	s.reuseAsNonZeroed(n)
 	var restore func()
 	if a == s {
 		s, restore = s.isolatedWorkspace(a)
@@ -578,7 +633,7 @@ func (s *SymDense) GrowSym(n int) Symmetric {
 // or the Eigendecomposition is not successful.
 func (s *SymDense) PowPSD(a Symmetric, pow float64) error {
 	dim := a.Symmetric()
-	s.reuseAs(dim)
+	s.reuseAsNonZeroed(dim)
 
 	var eigen EigenSym
 	ok := eigen.Factorize(a, true)

mat/triangular.go

@@ -15,6 +15,8 @@ import (
 var (
 	triDense *TriDense
 	_        Matrix            = triDense
+	_        allMatrix         = triDense
+	_        denseMatrix       = triDense
 	_        Triangular        = triDense
 	_        RawTriangular     = triDense
 	_        MutableTriangular = triDense
@@ -233,6 +235,7 @@ func (t *TriDense) SetRawTriangular(mat blas64.Triangular) {
 // Reset zeros the dimensions of the matrix so that it can be reused as the
 // receiver of a dimensionally restricted operation.
 //
+// Reset should not be used when the matrix shares backing data.
 // See the Reseter interface for more information.
 func (t *TriDense) Reset() {
 	// N and Stride must be zeroed in unison.
@@ -274,10 +277,32 @@ func untransposeTri(a Triangular) (Triangular, bool) {
 	return a, false
 }
 
-// reuseAs resizes a zero receiver to an n×n triangular matrix with the given
-// orientation. If the receiver is non-zero, reuseAs checks that the receiver
+// ReuseAsTri changes the receiver if it IsZero() to be of size n×n.
+//
+// ReuseAsTri re-uses the backing data slice if it has sufficient capacity,
+// otherwise a new slice is allocated. The data is then zeroed.
+//
+// ReuseAsTri panics if the receiver is not zero-sized, and panics if
+// the input size is less than one. To zero the receiver for re-use,
+// Reset should be used.
+func (t *TriDense) ReuseAsTri(n int, kind TriKind) {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
+		panic(ErrNegativeDimension)
+	}
+	if !t.IsZero() {
+		panic(ErrReuseNonZero)
+	}
+	t.reuseAsZeroed(n, kind)
+}
+
+// reuseAsNonZeroed resizes a zero receiver to an n×n triangular matrix with the given
+// orientation. If the receiver is non-zero, reuseAsNonZeroed checks that the receiver
 // is the correct size and orientation.
-func (t *TriDense) reuseAs(n int, kind TriKind) {
+func (t *TriDense) reuseAsNonZeroed(n int, kind TriKind) {
+	// reuseAsNonZeroed must be kept in sync with reuseAsZeroed.
 	if n == 0 {
 		panic(ErrZeroLength)
 	}
@@ -307,6 +332,41 @@ func (t *TriDense) reuseAs(n int, kind TriKind) {
 	}
 }
 
+// reuseAsZeroed resizes a zero receiver to an n×n triangular matrix with the given
+// orientation. If the receiver is non-zero, reuseAsZeroed checks that the receiver
+// is the correct size and orientation. It then zeros out the matrix data.
+func (t *TriDense) reuseAsZeroed(n int, kind TriKind) {
+	// reuseAsZeroed must be kept in sync with reuseAsNonZeroed.
+	if n == 0 {
+		panic(ErrZeroLength)
+	}
+	ul := blas.Lower
+	if kind == Upper {
+		ul = blas.Upper
+	}
+	if t.mat.N > t.cap {
+		panic(badTriCap)
+	}
+	if t.IsZero() {
+		t.mat = blas64.Triangular{
+			N:      n,
+			Stride: n,
+			Diag:   blas.NonUnit,
+			Data:   useZeroed(t.mat.Data, n*n),
+			Uplo:   ul,
+		}
+		t.cap = n
+		return
+	}
+	if t.mat.N != n {
+		panic(ErrShape)
+	}
+	if t.mat.Uplo != ul {
+		panic(ErrTriangle)
+	}
+	t.Zero()
+}
+
 // isolatedWorkspace returns a new TriDense matrix w with the size of a and
 // returns a callback to defer which performs cleanup at the return of the call.
 // This should be used when a method receiver is the same pointer as an input argument.
@@ -406,7 +466,7 @@ func (t *TriDense) Copy(a Matrix) (r, c int) {
 func (t *TriDense) InverseTri(a Triangular) error {
 	t.checkOverlapMatrix(a)
 	n, _ := a.Triangle()
-	t.reuseAs(a.Triangle())
+	t.reuseAsNonZeroed(a.Triangle())
 	t.Copy(a)
 	work := getFloats(3*n, false)
 	iwork := getInts(n, false)
@@ -443,7 +503,7 @@ func (t *TriDense) MulTri(a, b Triangular) {
 	bU, _ := untransposeTri(b)
 	t.checkOverlapMatrix(bU)
 	t.checkOverlapMatrix(aU)
-	t.reuseAs(n, kind)
+	t.reuseAsNonZeroed(n, kind)
 	var restore func()
 	if t == aU {
 		t, restore = t.isolatedWorkspace(aU)
@@ -506,7 +566,7 @@ func (t *TriDense) MulTri(a, b Triangular) {
 // the input, or ScaleTri will panic.
 func (t *TriDense) ScaleTri(f float64, a Triangular) {
 	n, kind := a.Triangle()
-	t.reuseAs(n, kind)
+	t.reuseAsNonZeroed(n, kind)
 
 	// TODO(btracey): Improve the set of fast-paths.
 	switch a := a.(type) {

mat/triband.go

@@ -16,6 +16,8 @@ var (
 
 	triBandDense *TriBandDense
 	_            Matrix           = triBandDense
+	_            allMatrix        = triBandDense
+	_            denseMatrix      = triBandDense
 	_            Triangular       = triBandDense
 	_            Banded           = triBandDense
 	_            TriBanded        = triBandDense

mat/vector.go

@@ -13,9 +13,10 @@ import (
 var (
 	vector *VecDense
 
-	_ Matrix  = vector
-	_ Vector  = vector
-	_ Reseter = vector
+	_ Matrix    = vector
+	_ allMatrix = vector
+	_ Vector    = vector
+	_ Reseter   = vector
 )
 
 // Vector is a vector.
@@ -252,7 +253,7 @@ func (v *VecDense) ScaleVec(alpha float64, a Vector) {
 		return
 	}
 
-	v.reuseAs(n)
+	v.reuseAsNonZeroed(n)
 
 	if rv, ok := a.(RawVectorer); ok {
 		mat := rv.RawVector()
@@ -310,7 +311,7 @@ func (v *VecDense) AddScaledVec(a Vector, alpha float64, b Vector) {
 		fast = false
 	}
 
-	v.reuseAs(ar)
+	v.reuseAsNonZeroed(ar)
 
 	switch {
 	case alpha == 0: // v <- a
@@ -353,7 +354,7 @@ func (v *VecDense) AddVec(a, b Vector) {
 		panic(ErrShape)
 	}
 
-	v.reuseAs(ar)
+	v.reuseAsNonZeroed(ar)
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
@@ -396,7 +397,7 @@ func (v *VecDense) SubVec(a, b Vector) {
 		panic(ErrShape)
 	}
 
-	v.reuseAs(ar)
+	v.reuseAsNonZeroed(ar)
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
@@ -440,7 +441,7 @@ func (v *VecDense) MulElemVec(a, b Vector) {
 		panic(ErrShape)
 	}
 
-	v.reuseAs(ar)
+	v.reuseAsNonZeroed(ar)
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
@@ -489,7 +490,7 @@ func (v *VecDense) DivElemVec(a, b Vector) {
 		panic(ErrShape)
 	}
 
-	v.reuseAs(ar)
+	v.reuseAsNonZeroed(ar)
 
 	aU, _ := untransposeExtract(a)
 	bU, _ := untransposeExtract(b)
@@ -550,7 +551,7 @@ func (v *VecDense) MulVec(a Matrix, b Vector) {
 		fast = false
 	}
 
-	v.reuseAs(r)
+	v.reuseAsNonZeroed(r)
 	var restore func()
 	if v == aU {
 		v, restore = v.isolatedWorkspace(aU.(*VecDense))
@@ -639,9 +640,31 @@ func (v *VecDense) MulVec(a Matrix, b Vector) {
 	}
 }
 
-// reuseAs resizes an empty vector to a r×1 vector,
+// ReuseAsVec changes the receiver if it IsZero() to be of size n×1.
+//
+// ReuseAsVec re-uses the backing data slice if it has sufficient capacity,
+// otherwise a new slice is allocated. The data is then zeroed.
+//
+// ReuseAsVec panics if the receiver is not zero-sized, and panics if
+// the input size is less than one. To zero the receiver for re-use,
+// Reset should be used.
+func (v *VecDense) ReuseAsVec(n int) {
+	if n <= 0 {
+		if n == 0 {
+			panic(ErrZeroLength)
+		}
+		panic(ErrNegativeDimension)
+	}
+	if !v.IsZero() {
+		panic(ErrReuseNonZero)
+	}
+	v.reuseAsZeroed(n)
+}
+
+// reuseAsNonZeroed resizes an empty vector to a r×1 vector,
 // or checks that a non-empty matrix is r×1.
-func (v *VecDense) reuseAs(r int) {
+func (v *VecDense) reuseAsNonZeroed(r int) {
+	// reuseAsNonZeroed must be kept in sync with reuseAsZeroed.
 	if r == 0 {
 		panic(ErrZeroLength)
 	}
@@ -658,6 +681,27 @@ func (v *VecDense) reuseAs(r int) {
 	}
 }
 
+// reuseAsZeroed resizes an empty vector to a r×1 vector,
+// or checks that a non-empty matrix is r×1.
+func (v *VecDense) reuseAsZeroed(r int) {
+	// reuseAsZeroed must be kept in sync with reuseAsNonZeroed.
+	if r == 0 {
+		panic(ErrZeroLength)
+	}
+	if v.IsZero() {
+		v.mat = blas64.Vector{
+			N:    r,
+			Inc:  1,
+			Data: useZeroed(v.mat.Data, r),
+		}
+		return
+	}
+	if r != v.mat.N {
+		panic(ErrShape)
+	}
+	v.Zero()
+}
+
 // IsZero returns whether the receiver is zero-sized. Zero-sized vectors can be the
 // receiver for size-restricted operations. VecDenses can be zeroed using Reset.
 func (v *VecDense) IsZero() bool {