mat: Add Zero method to reset values of matrices (#819)

* mat: Add Zero method to reset values of matrices

mat/band.go

@@ -242,3 +242,15 @@ func (b *BandDense) DoColNonZero(j int, fn func(i, j int, v float64)) {
 		}
 	}
 }
+
+// Zero sets all of the matrix elements to zero.
+func (b *BandDense) Zero() {
+	m := b.mat.Rows
+	kL := b.mat.KL
+	nCol := b.mat.KU + 1 + kL
+	for i := 0; i < m; i++ {
+		l := max(0, kL-i)
+		u := min(nCol, m+kL-i)
+		zero(b.mat.Data[i*b.mat.Stride+l : i*b.mat.Stride+u])
+	}
+}

mat/band_test.go

@@ -248,6 +248,59 @@ func TestNewDiagonalRect(t *testing.T) {
 	}
 }
 
+func TestBandDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*BandDense{
+		&BandDense{
+			mat: blas64.Band{
+				Rows:   6,
+				Cols:   7,
+				Stride: 8,
+				KL:     1,
+				KU:     2,
+				Data: []float64{
+					-1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, -1, -1, -1, -1, -1,
+					1, 1, -1, -1, -1, -1, -1, -1,
+				},
+			},
+		},
+		&BandDense{
+			mat: blas64.Band{
+				Rows:   6,
+				Cols:   7,
+				Stride: 8,
+				KL:     2,
+				KU:     1,
+				Data: []float64{
+					-1, -1, 1, 1, -1, -1, -1, -1,
+					-1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, 1, -1, -1, -1, -1,
+					1, 1, 1, -1, -1, -1, -1, -1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}
+
 func TestBandDiagView(t *testing.T) {
 	for cas, test := range []*BandDense{
 		NewBandDense(1, 1, 0, 0, []float64{1}),

mat/dense.go

@@ -122,6 +122,13 @@ func (m *Dense) reuseAsZeroed(r, c int) {
 	if r != m.mat.Rows || c != m.mat.Cols {
 		panic(ErrShape)
 	}
+	m.Zero()
+}
+
+// Zero sets all of the matrix elements to zero.
+func (m *Dense) Zero() {
+	r := m.mat.Rows
+	c := m.mat.Cols
 	for i := 0; i < r; i++ {
 		zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c])
 	}

mat/dense_test.go

@@ -257,6 +257,38 @@ func TestSetRowColumn(t *testing.T) {
 	}
 }
 
+func TestDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*Dense{
+		&Dense{
+			mat: blas64.General{
+				Rows:   4,
+				Cols:   3,
+				Stride: 5,
+				Data: []float64{
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}
+
 func TestRowColView(t *testing.T) {
 	for _, test := range []struct {
 		mat [][]float64

mat/diagonal.go

@@ -168,6 +168,13 @@ func (d *DiagDense) Reset() {
 	d.mat.Data = d.mat.Data[:0]
 }
 
+// Zero sets all of the matrix elements to zero.
+func (d *DiagDense) Zero() {
+	for i := 0; i < d.mat.N; i++ {
+		d.mat.Data[d.mat.Inc*i] = 0
+	}
+}
+
 // DiagView returns the diagonal as a matrix backed by the original data.
 func (d *DiagDense) DiagView() Diagonal {
 	return d

mat/diagonal_test.go

@@ -56,6 +56,38 @@ func TestNewDiagDense(t *testing.T) {
 	}
 }
 
+func TestDiagDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*DiagDense{
+		&DiagDense{
+			mat: blas64.Vector{
+				N:   5,
+				Inc: 2,
+				Data: []float64{
+					1, -1,
+					1, -1,
+					1, -1,
+					1, -1,
+					1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}
+
 func TestDiagonalStride(t *testing.T) {
 	for _, test := range []struct {
 		diag  *DiagDense

mat/symband.go

@@ -145,6 +145,14 @@ func (s *SymBandDense) RawSymBand() blas64.SymmetricBand {
 	return s.mat
 }
 
+// Zero sets all of the matrix elements to zero.
+func (s *SymBandDense) Zero() {
+	for i := 0; i < s.mat.N; i++ {
+		u := min(1+s.mat.K, s.mat.N-i)
+		zero(s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+u])
+	}
+}
+
 // DiagView returns the diagonal as a matrix backed by the original data.
 func (s *SymBandDense) DiagView() Diagonal {
 	n := s.mat.N

mat/symband_test.go

@@ -195,3 +195,38 @@ func TestSymBandDiagView(t *testing.T) {
 		testDiagView(t, cas, test)
 	}
 }
+
+func TestSymBandDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*SymBandDense{
+		&SymBandDense{
+			mat: blas64.SymmetricBand{
+				Uplo:   blas.Upper,
+				N:      6,
+				K:      2,
+				Stride: 5,
+				Data: []float64{
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, -1, -1, -1,
+					1, -1, -1, -1, -1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}

mat/symmetric.go

@@ -131,6 +131,13 @@ func (s *SymDense) Reset() {
 	s.mat.Data = s.mat.Data[:0]
 }
 
+// Zero sets all of the matrix elements to zero.
+func (s *SymDense) Zero() {
+	for i := 0; i < s.mat.N; i++ {
+		zero(s.mat.Data[i*s.mat.Stride+i : i*s.mat.Stride+s.mat.N])
+	}
+}
+
 // IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
 // receiver for size-restricted operations. SymDense matrices can be zeroed using Reset.
 func (s *SymDense) IsZero() bool {

mat/symmetric_test.go

@@ -129,6 +129,38 @@ func TestSymAtSet(t *testing.T) {
 	}
 }
 
+func TestSymDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*SymDense{
+		&SymDense{
+			mat: blas64.Symmetric{
+				Uplo:   blas.Upper,
+				N:      4,
+				Stride: 5,
+				Data: []float64{
+					1, 1, 1, 1, -1,
+					-1, 1, 1, 1, -1,
+					-1, -1, 1, 1, -1,
+					-1, -1, -1, 1, -1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}
+
 func TestSymDiagView(t *testing.T) {
 	for cas, test := range []*SymDense{
 		NewSymDense(1, []float64{1}),

mat/triangular.go

@@ -230,6 +230,19 @@ func (t *TriDense) Reset() {
 	t.mat.Data = t.mat.Data[:0]
 }
 
+// Zero sets all of the matrix elements to zero.
+func (t *TriDense) Zero() {
+	if t.isUpper() {
+		for i := 0; i < t.mat.N; i++ {
+			zero(t.mat.Data[i*t.mat.Stride+i : i*t.mat.Stride+t.mat.N])
+		}
+		return
+	}
+	for i := 0; i < t.mat.N; i++ {
+		zero(t.mat.Data[i*t.mat.Stride : i*t.mat.Stride+i+1])
+	}
+}
+
 // IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
 // receiver for size-restricted operations. TriDense matrices can be zeroed using Reset.
 func (t *TriDense) IsZero() bool {

mat/triangular_test.go

@@ -139,6 +139,51 @@ func TestTriAtSet(t *testing.T) {
 	}
 }
 
+func TestTriDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*TriDense{
+		&TriDense{
+			mat: blas64.Triangular{
+				Uplo:   blas.Upper,
+				N:      4,
+				Stride: 5,
+				Data: []float64{
+					1, 1, 1, 1, -1,
+					-1, 1, 1, 1, -1,
+					-1, -1, 1, 1, -1,
+					-1, -1, -1, 1, -1,
+				},
+			},
+		},
+		&TriDense{
+			mat: blas64.Triangular{
+				Uplo:   blas.Lower,
+				N:      4,
+				Stride: 5,
+				Data: []float64{
+					1, -1, -1, -1, -1,
+					1, 1, -1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, 1, -1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}
+
 func TestTriDiagView(t *testing.T) {
 	for cas, test := range []*TriDense{
 		NewTriDense(1, Upper, []float64{1}),

mat/triband.go

@@ -256,6 +256,21 @@ func (t *TriBandDense) Reset() {
 	t.mat.Data = t.mat.Data[:0]
 }
 
+// Zero sets all of the matrix elements to zero.
+func (t *TriBandDense) Zero() {
+	if t.isUpper() {
+		for i := 0; i < t.mat.N; i++ {
+			u := min(1+t.mat.K, t.mat.N-i)
+			zero(t.mat.Data[i*t.mat.Stride : i*t.mat.Stride+u])
+		}
+		return
+	}
+	for i := 0; i < t.mat.N; i++ {
+		l := max(0, t.mat.K-i)
+		zero(t.mat.Data[i*t.mat.Stride+l : i*t.mat.Stride+t.mat.K+1])
+	}
+}
+
 func (t *TriBandDense) isUpper() bool {
 	return isUpperUplo(t.mat.Uplo)
 }

mat/triband_test.go

@@ -333,6 +333,57 @@ func TestTriBandAtSetUpper(t *testing.T) {
 	}
 }
 
+func TestTriBandDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*TriBandDense{
+		&TriBandDense{
+			mat: blas64.TriangularBand{
+				Uplo:   blas.Upper,
+				N:      6,
+				K:      2,
+				Stride: 5,
+				Data: []float64{
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, -1, -1, -1,
+					1, -1, -1, -1, -1,
+				},
+			},
+		},
+		&TriBandDense{
+			mat: blas64.TriangularBand{
+				Uplo:   blas.Lower,
+				N:      6,
+				K:      2,
+				Stride: 5,
+				Data: []float64{
+					-1, -1, 1, -1, -1,
+					-1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+					1, 1, 1, -1, -1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}
+
 func TestTriBandDiagView(t *testing.T) {
 	for cas, test := range []*TriBandDense{
 		NewTriBandDense(1, 0, Upper, []float64{1}),

mat/vector.go

@@ -178,6 +178,13 @@ func (v *VecDense) Reset() {
 	v.mat.Data = v.mat.Data[:0]
 }
 
+// Zero sets all of the matrix elements to zero.
+func (v *VecDense) Zero() {
+	for i := 0; i < v.mat.N; i++ {
+		v.mat.Data[v.mat.Inc*i] = 0
+	}
+}
+
 // CloneVec makes a copy of a into the receiver, overwriting the previous value
 // of the receiver.
 func (v *VecDense) CloneVec(a Vector) {

mat/vector_test.go

@@ -181,6 +181,38 @@ func TestVecDenseAtSet(t *testing.T) {
 	}
 }
 
+func TestVecDenseZero(t *testing.T) {
+	// Elements that equal 1 should be set to zero, elements that equal -1
+	// should remain unchanged.
+	for _, test := range []*VecDense{
+		&VecDense{
+			mat: blas64.Vector{
+				N:   5,
+				Inc: 2,
+				Data: []float64{
+					1, -1,
+					1, -1,
+					1, -1,
+					1, -1,
+					1,
+				},
+			},
+		},
+	} {
+		dataCopy := make([]float64, len(test.mat.Data))
+		copy(dataCopy, test.mat.Data)
+		test.Zero()
+		for i, v := range test.mat.Data {
+			if dataCopy[i] != -1 && v != 0 {
+				t.Errorf("Matrix not zeroed in bounds")
+			}
+			if dataCopy[i] == -1 && v != -1 {
+				t.Errorf("Matrix zeroed out of bounds")
+			}
+		}
+	}
+}
+
 func TestVecDenseMul(t *testing.T) {
 	method := func(receiver, a, b Matrix) {
 		type mulVecer interface {