Merge 44689b6 into d0fc9ef

matrix/mat64/gsvd.go

@@ -203,9 +203,10 @@ func (gsvd *GSVD) ValuesB(s []float64) []float64 {
 	return s
 }
 
-// ZeroRFromGSVD extracts the matrix [ 0 R ] from the singular value decomposition, storing
-// the result in-place into the receiver. [ 0 R ] is size (k+l)×c.
-func (m *Dense) ZeroRFromGSVD(gsvd *GSVD) {
+// ZeroRTo extracts the matrix [ 0 R ] from the singular value decomposition, storing
+// the result in-place into dst. [ 0 R ] is size (k+l)×c.
+// If dst is nil, a new matrix is allocated. The resulting ZeroR matrix is returned.
+func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense {
 	if gsvd.kind == 0 {
 		panic("gsvd: no decomposition computed")
 	}
@@ -214,112 +215,147 @@ func (m *Dense) ZeroRFromGSVD(gsvd *GSVD) {
 	k := gsvd.k
 	l := gsvd.l
 	h := min(k+l, r)
-	m.reuseAsZeroed(k+l, c)
+	if dst == nil {
+		dst = NewDense(k+l, c, nil)
+	} else {
+		dst.reuseAsZeroed(k+l, c)
+	}
 	a := Dense{
 		mat:     gsvd.a,
 		capRows: r,
 		capCols: c,
 	}
-	m.Slice(0, h, c-k-l, c).(*Dense).
+	dst.Slice(0, h, c-k-l, c).(*Dense).
 		Copy(a.Slice(0, h, c-k-l, c))
 	if r < k+l {
 		b := Dense{
 			mat:     gsvd.b,
 			capRows: gsvd.p,
 			capCols: c,
 		}
-		m.Slice(r, k+l, c+r-k-l, c).(*Dense).
+		dst.Slice(r, k+l, c+r-k-l, c).(*Dense).
 			Copy(b.Slice(r-k, l, c+r-k-l, c))
 	}
+	return dst
 }
 
-// SigmaAFromGSVD extracts the matrix Σ₁ from the singular value decomposition, storing
-// the result in-place into the receiver. Σ₁ is size r×(k+l).
-func (m *Dense) SigmaAFromGSVD(gsvd *GSVD) {
+// SigmaATo extracts the matrix Σ₁ from the singular value decomposition, storing
+// the result in-place into dst. Σ₁ is size r×(k+l).
+// If dst is nil, a new matrix is allocated. The resulting SigmaA matrix is returned.
+func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense {
 	if gsvd.kind == 0 {
 		panic("gsvd: no decomposition computed")
 	}
 	r := gsvd.r
 	k := gsvd.k
 	l := gsvd.l
-	m.reuseAsZeroed(r, k+l)
+	if dst == nil {
+		dst = NewDense(r, k+l, nil)
+	} else {
+		dst.reuseAsZeroed(r, k+l)
+	}
 	for i := 0; i < k; i++ {
-		m.set(i, i, 1)
+		dst.set(i, i, 1)
 	}
 	for i := k; i < min(r, k+l); i++ {
-		m.set(i, i, gsvd.s1[i])
+		dst.set(i, i, gsvd.s1[i])
 	}
+	return dst
 }
 
-// SigmaBFromGSVD extracts the matrix Σ₂ from the singular value decomposition, storing
-// the result in-place into the receiver. Σ₂ is size p×(k+l).
-func (m *Dense) SigmaBFromGSVD(gsvd *GSVD) {
+// SigmaBTo extracts the matrix Σ₂ from the singular value decomposition, storing
+// the result in-place into dst. Σ₂ is size p×(k+l).
+// If dst is nil, a new matrix is allocated. The resulting SigmaB matrix is returned.
+func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense {
 	if gsvd.kind == 0 {
 		panic("gsvd: no decomposition computed")
 	}
 	r := gsvd.r
 	p := gsvd.p
 	k := gsvd.k
 	l := gsvd.l
-	m.reuseAsZeroed(p, k+l)
+	if dst == nil {
+		dst = NewDense(p, k+l, nil)
+	} else {
+		dst.reuseAsZeroed(p, k+l)
+	}
 	for i := 0; i < min(l, r-k); i++ {
-		m.set(i, i+k, gsvd.s2[k+i])
+		dst.set(i, i+k, gsvd.s2[k+i])
 	}
 	for i := r - k; i < l; i++ {
-		m.set(i, i+k, 1)
+		dst.set(i, i+k, 1)
 	}
+	return dst
 }
 
-// UFromGSVD extracts the matrix U from the singular value decomposition, storing
-// the result in-place into the receiver. U is size r×r.
-func (m *Dense) UFromGSVD(gsvd *GSVD) {
+// UTo extracts the matrix U from the singular value decomposition, storing
+// the result in-place into dst. U is size r×r.
+// If dst is nil, a new matrix is allocated. The resulting U matrix is returned.
+func (gsvd *GSVD) UTo(dst *Dense) *Dense {
 	if gsvd.kind&matrix.GSVDU == 0 {
 		panic("mat64: improper GSVD kind")
 	}
 	r := gsvd.u.Rows
 	c := gsvd.u.Cols
-	m.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     gsvd.u,
 		capRows: r,
 		capCols: c,
 	}
-	m.Copy(tmp)
+	dst.Copy(tmp)
+	return dst
 }
 
-// VFromGSVD extracts the matrix V from the singular value decomposition, storing
-// the result in-place into the receiver. V is size p×p.
-func (m *Dense) VFromGSVD(gsvd *GSVD) {
+// VTo extracts the matrix V from the singular value decomposition, storing
+// the result in-place into dst. V is size p×p.
+// If dst is nil, a new matrix is allocated. The resulting V matrix is returned.
+func (gsvd *GSVD) VTo(dst *Dense) *Dense {
 	if gsvd.kind&matrix.GSVDV == 0 {
 		panic("mat64: improper GSVD kind")
 	}
 	r := gsvd.v.Rows
 	c := gsvd.v.Cols
-	m.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     gsvd.v,
 		capRows: r,
 		capCols: c,
 	}
-	m.Copy(tmp)
+	dst.Copy(tmp)
+	return dst
 }
 
-// QFromGSVD extracts the matrix Q from the singular value decomposition, storing
-// the result in-place into the receiver. Q is size c×c.
-func (m *Dense) QFromGSVD(gsvd *GSVD) {
+// QTo extracts the matrix Q from the singular value decomposition, storing
+// the result in-place into dst. Q is size c×c.
+// If dst is nil, a new matrix is allocated. The resulting Q matrix is returned.
+func (gsvd *GSVD) QTo(dst *Dense) *Dense {
 	if gsvd.kind&matrix.GSVDQ == 0 {
 		panic("mat64: improper GSVD kind")
 	}
 	r := gsvd.q.Rows
 	c := gsvd.q.Cols
-	m.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     gsvd.q,
 		capRows: r,
 		capCols: c,
 	}
-	m.Copy(tmp)
+	dst.Copy(tmp)
+	return dst
 }

matrix/mat64/gsvd_example_test.go

@@ -25,22 +25,19 @@ func ExampleGSVD() {
 		log.Fatal("GSVD factorization failed")
 	}
 
-	var u, v mat64.Dense
-	u.UFromGSVD(&gsvd)
-	v.VFromGSVD(&gsvd)
+	u := gsvd.UTo(nil)
+	v := gsvd.VTo(nil)
 
 	s1 := gsvd.ValuesA(nil)
 	s2 := gsvd.ValuesB(nil)
 
 	fmt.Printf("Africa\n\ts1 = %.4f\n\n\tU = %.4f\n\n",
-		s1, mat64.Formatted(&u, mat64.Prefix("\t    "), mat64.Excerpt(2)))
+		s1, mat64.Formatted(u, mat64.Prefix("\t    "), mat64.Excerpt(2)))
 	fmt.Printf("Latin America/Caribbean\n\ts2 = %.4f\n\n\tV = %.4f\n",
-		s2, mat64.Formatted(&v, mat64.Prefix("\t    "), mat64.Excerpt(2)))
+		s2, mat64.Formatted(v, mat64.Prefix("\t    "), mat64.Excerpt(2)))
 
-	var zeroR, q mat64.Dense
-	zeroR.ZeroRFromGSVD(&gsvd)
-	q.QFromGSVD(&gsvd)
-	q.Mul(&zeroR, &q)
+	var q mat64.Dense
+	q.Mul(gsvd.ZeroRTo(nil), gsvd.QTo(nil))
 	fmt.Printf("\nCommon basis vectors\n\n\tQ^T = %.4f\n",
 		mat64.Formatted(q.T(), mat64.Prefix("\t      ")))
 

matrix/mat64/gsvd_test.go

@@ -106,14 +106,13 @@ func TestGSVD(t *testing.T) {
 }
 
 func extractGSVD(gsvd *GSVD) (c, s []float64, s1, s2, zR, u, v, q *Dense) {
-	var s1m, s2m, zeroR, um, vm, qm Dense
-	s1m.SigmaAFromGSVD(gsvd)
-	s2m.SigmaBFromGSVD(gsvd)
-	zeroR.ZeroRFromGSVD(gsvd)
-	um.UFromGSVD(gsvd)
-	vm.VFromGSVD(gsvd)
-	qm.QFromGSVD(gsvd)
+	s1 = gsvd.SigmaATo(nil)
+	s2 = gsvd.SigmaBTo(nil)
+	zR = gsvd.ZeroRTo(nil)
+	u = gsvd.UTo(nil)
+	v = gsvd.VTo(nil)
+	q = gsvd.QTo(nil)
 	c = gsvd.ValuesA(nil)
 	s = gsvd.ValuesB(nil)
-	return c, s, &s1m, &s2m, &zeroR, &um, &vm, &qm
+	return c, s, s1, s2, zR, u, v, q
 }

matrix/mat64/hogsvd.go

@@ -141,24 +141,31 @@ func (gsvd *HOGSVD) Len() int {
 	return gsvd.n
 }
 
-// UFromHOGSVD extracts the matrix U_n from the singular value decomposition, storing
-// the result in-place into the receiver. U_n is size r×c.
+// UTo extracts the matrix U_n from the singular value decomposition, storing
+// the result in-place into dst. U_n is size r×c.
+// If dst is nil, a new matrix is allocated. The resulting U matrix is returned.
 //
-// UFromHOGSVD will panic if the receiver does not contain a successful factorization.
-func (m *Dense) UFromHOGSVD(gsvd *HOGSVD, n int) {
+// UTo will panic if the receiver does not contain a successful factorization.
+func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense {
 	if gsvd.n == 0 {
 		panic("hogsvd: unsuccessful factorization")
 	}
 	if n < 0 || gsvd.n <= n {
 		panic("hogsvd: invalid index")
 	}
 
-	m.reuseAs(gsvd.b[n].Dims())
-	m.Copy(&gsvd.b[n])
+	if dst == nil {
+		r, c := gsvd.b[n].Dims()
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(gsvd.b[n].Dims())
+	}
+	dst.Copy(&gsvd.b[n])
 	for j, f := range gsvd.Values(nil, n) {
-		v := m.ColView(j)
+		v := dst.ColView(j)
 		v.ScaleVec(1/f, v)
 	}
+	return dst
 }
 
 // Values returns the nth set of singular values of the factorized system.
@@ -188,13 +195,21 @@ func (gsvd *HOGSVD) Values(s []float64, n int) []float64 {
 	return s
 }
 
-// VFromHOGSVD extracts the matrix V from the singular value decomposition, storing
-// the result in-place into the receiver. V is size c×c.
+// VTo extracts the matrix V from the singular value decomposition, storing
+// the result in-place into dst. V is size c×c.
+// If dst is nil, a new matrix is allocated. The resulting V matrix is returned.
 //
-// VFromHOGSVD will panic if the receiver does not contain a successful factorization.
-func (m *Dense) VFromHOGSVD(gsvd *HOGSVD) {
+// VTo will panic if the receiver does not contain a successful factorization.
+func (gsvd *HOGSVD) VTo(dst *Dense) *Dense {
 	if gsvd.n == 0 {
 		panic("hogsvd: unsuccessful factorization")
 	}
-	*m = *DenseCopyOf(gsvd.v)
+	if dst == nil {
+		r, c := gsvd.v.Dims()
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(gsvd.v.Dims())
+	}
+	dst.Copy(gsvd.v)
+	return dst
 }

matrix/mat64/hogsvd_example_test.go

@@ -24,15 +24,13 @@ func ExampleHOGSVD() {
 	}
 
 	for i, n := range []string{"Africa", "Asia", "Latin America/Caribbean", "Oceania"} {
-		var u mat64.Dense
-		u.UFromHOGSVD(&gsvd, i)
+		u := gsvd.UTo(nil, i)
 		s := gsvd.Values(nil, i)
 		fmt.Printf("%s\n\ts_%d = %.4f\n\n\tU_%[2]d = %.4[4]f\n",
-			n, i, s, mat64.Formatted(&u, mat64.Prefix("\t      ")))
+			n, i, s, mat64.Formatted(u, mat64.Prefix("\t      ")))
 	}
 
-	var v mat64.Dense
-	v.VFromHOGSVD(&gsvd)
+	v := gsvd.VTo(nil)
 	fmt.Printf("\nCommon basis vectors\n\n\tV^T = %.4f",
 		mat64.Formatted(v.T(), mat64.Prefix("\t      ")))
 

matrix/mat64/hogsvd_test.go

@@ -79,11 +79,9 @@ func extractHOGSVD(gsvd *HOGSVD) (u []*Dense, s [][]float64, v *Dense) {
 	u = make([]*Dense, gsvd.Len())
 	s = make([][]float64, gsvd.Len())
 	for i := 0; i < gsvd.Len(); i++ {
-		u[i] = &Dense{}
-		u[i].UFromHOGSVD(gsvd, i)
+		u[i] = gsvd.UTo(nil, i)
 		s[i] = gsvd.Values(nil, i)
 	}
-	v = &Dense{}
-	v.VFromHOGSVD(gsvd)
+	v = gsvd.VTo(nil)
 	return u, s, v
 }