diff --git a/matrix/mat64/gsvd.go b/matrix/mat64/gsvd.go
index e3064848..a87a94d6 100644
--- a/matrix/mat64/gsvd.go
+++ b/matrix/mat64/gsvd.go
@@ -205,7 +205,8 @@ func (gsvd *GSVD) ValuesB(s []float64) []float64 {
 
 // 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.
-func (gsvd *GSVD) ZeroRTo(dst *Dense) {
+// 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,7 +215,11 @@ func (gsvd *GSVD) ZeroRTo(dst *Dense) {
 	k := gsvd.k
 	l := gsvd.l
 	h := min(k+l, r)
-	dst.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,
@@ -231,29 +236,37 @@ func (gsvd *GSVD) ZeroRTo(dst *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
 }
 
 // SigmaATo extracts the matrix Σ₁ from the singular value decomposition, storing
 // the result in-place into dst. Σ₁ is size r×(k+l).
-func (gsvd *GSVD) SigmaATo(dst *Dense) {
+// 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
-	dst.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++ {
 		dst.set(i, i, 1)
 	}
 	for i := k; i < min(r, k+l); i++ {
 		dst.set(i, i, gsvd.s1[i])
 	}
+	return dst
 }
 
 // SigmaBTo extracts the matrix Σ₂ from the singular value decomposition, storing
 // the result in-place into dst. Σ₂ is size p×(k+l).
-func (gsvd *GSVD) SigmaBTo(dst *Dense) {
+// 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")
 	}
@@ -261,24 +274,34 @@ func (gsvd *GSVD) SigmaBTo(dst *Dense) {
 	p := gsvd.p
 	k := gsvd.k
 	l := gsvd.l
-	dst.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++ {
 		dst.set(i, i+k, gsvd.s2[k+i])
 	}
 	for i := r - k; i < l; i++ {
 		dst.set(i, i+k, 1)
 	}
+	return dst
 }
 
 // UTo extracts the matrix U from the singular value decomposition, storing
 // the result in-place into dst. U is size r×r.
-func (gsvd *GSVD) UTo(dst *Dense) {
+// 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
-	dst.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     gsvd.u,
@@ -286,17 +309,23 @@ func (gsvd *GSVD) UTo(dst *Dense) {
 		capCols: c,
 	}
 	dst.Copy(tmp)
+	return dst
 }
 
 // VTo extracts the matrix V from the singular value decomposition, storing
 // the result in-place into dst. V is size p×p.
-func (gsvd *GSVD) VTo(dst *Dense) {
+// 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
-	dst.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     gsvd.v,
@@ -304,17 +333,23 @@ func (gsvd *GSVD) VTo(dst *Dense) {
 		capCols: c,
 	}
 	dst.Copy(tmp)
+	return dst
 }
 
 // QTo extracts the matrix Q from the singular value decomposition, storing
 // the result in-place into dst. Q is size c×c.
-func (gsvd *GSVD) QTo(dst *Dense) {
+// 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
-	dst.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     gsvd.q,
@@ -322,4 +357,5 @@ func (gsvd *GSVD) QTo(dst *Dense) {
 		capCols: c,
 	}
 	dst.Copy(tmp)
+	return dst
 }
diff --git a/matrix/mat64/gsvd_example_test.go b/matrix/mat64/gsvd_example_test.go
index 5f605642..6d34f712 100644
--- a/matrix/mat64/gsvd_example_test.go
+++ b/matrix/mat64/gsvd_example_test.go
@@ -25,22 +25,19 @@ func ExampleGSVD() {
 		log.Fatal("GSVD factorization failed")
 	}
 
-	var u, v mat64.Dense
-	gsvd.UTo(&u)
-	gsvd.VTo(&v)
+	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
-	gsvd.ZeroRTo(&zeroR)
-	gsvd.QTo(&q)
-	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      ")))
 
diff --git a/matrix/mat64/gsvd_test.go b/matrix/mat64/gsvd_test.go
index 0e4d8274..e8b82083 100644
--- a/matrix/mat64/gsvd_test.go
+++ b/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
-	gsvd.SigmaATo(&s1m)
-	gsvd.SigmaBTo(&s2m)
-	gsvd.ZeroRTo(&zeroR)
-	gsvd.UTo(&um)
-	gsvd.VTo(&vm)
-	gsvd.QTo(&qm)
+	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
 }
diff --git a/matrix/mat64/hogsvd.go b/matrix/mat64/hogsvd.go
index d9ea0b5d..9cd6ec4f 100644
--- a/matrix/mat64/hogsvd.go
+++ b/matrix/mat64/hogsvd.go
@@ -143,9 +143,10 @@ func (gsvd *HOGSVD) Len() int {
 
 // 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.
 //
 // UTo will panic if the receiver does not contain a successful factorization.
-func (gsvd *HOGSVD) UTo(dst *Dense, n int) {
+func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense {
 	if gsvd.n == 0 {
 		panic("hogsvd: unsuccessful factorization")
 	}
@@ -153,12 +154,18 @@ func (gsvd *HOGSVD) UTo(dst *Dense, n int) {
 		panic("hogsvd: invalid index")
 	}
 
-	dst.reuseAs(gsvd.b[n].Dims())
+	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 := dst.ColView(j)
 		v.ScaleVec(1/f, v)
 	}
+	return dst
 }
 
 // Values returns the nth set of singular values of the factorized system.
@@ -190,12 +197,19 @@ func (gsvd *HOGSVD) Values(s []float64, n int) []float64 {
 
 // 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.
 //
 // VTo will panic if the receiver does not contain a successful factorization.
-func (gsvd *HOGSVD) VTo(dst *Dense) {
+func (gsvd *HOGSVD) VTo(dst *Dense) *Dense {
 	if gsvd.n == 0 {
 		panic("hogsvd: unsuccessful factorization")
 	}
-	dst.reuseAs(gsvd.v.Dims())
+	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
 }
diff --git a/matrix/mat64/hogsvd_example_test.go b/matrix/mat64/hogsvd_example_test.go
index 8fb2f00a..dfb0821f 100644
--- a/matrix/mat64/hogsvd_example_test.go
+++ b/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
-		gsvd.UTo(&u, 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
-	gsvd.VTo(&v)
+	v := gsvd.VTo(nil)
 	fmt.Printf("\nCommon basis vectors\n\n\tV^T = %.4f",
 		mat64.Formatted(v.T(), mat64.Prefix("\t      ")))
 
diff --git a/matrix/mat64/hogsvd_test.go b/matrix/mat64/hogsvd_test.go
index 320fa2bc..f6543944 100644
--- a/matrix/mat64/hogsvd_test.go
+++ b/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{}
-		gsvd.UTo(u[i], i)
+		u[i] = gsvd.UTo(nil, i)
 		s[i] = gsvd.Values(nil, i)
 	}
-	v = &Dense{}
-	gsvd.VTo(v)
+	v = gsvd.VTo(nil)
 	return u, s, v
 }
diff --git a/matrix/mat64/lq.go b/matrix/mat64/lq.go
index f8a04732..b9bffd36 100644
--- a/matrix/mat64/lq.go
+++ b/matrix/mat64/lq.go
@@ -59,11 +59,16 @@ func (lq *LQ) Factorize(a Matrix) {
 // and upper triangular matrices.
 
 // LTo extracts the m×n lower trapezoidal matrix from a LQ decomposition.
-func (lq *LQ) LTo(dst *Dense) {
+// If dst is nil, a new matrix is allocated. The resulting L matrix is returned.
+func (lq *LQ) LTo(dst *Dense) *Dense {
 	r, c := lq.lq.Dims()
-	dst.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
-	// Disguise the LQ as a lower triangular
+	// Disguise the LQ as a lower triangular.
 	t := &TriDense{
 		mat: blas64.Triangular{
 			N:      r,
@@ -77,18 +82,25 @@ func (lq *LQ) LTo(dst *Dense) {
 	dst.Copy(t)
 
 	if r == c {
-		return
+		return dst
 	}
 	// Zero right of the triangular.
 	for i := 0; i < r; i++ {
 		zero(dst.mat.Data[i*dst.mat.Stride+r : i*dst.mat.Stride+c])
 	}
+
+	return dst
 }
 
 // QTo extracts the n×n orthonormal matrix Q from an LQ decomposition.
-func (lq *LQ) QTo(dst *Dense) {
+// If dst is nil, a new matrix is allocated. The resulting Q matrix is returned.
+func (lq *LQ) QTo(dst *Dense) *Dense {
 	r, c := lq.lq.Dims()
-	dst.reuseAs(c, c)
+	if dst == nil {
+		dst = NewDense(c, c, nil)
+	} else {
+		dst.reuseAs(c, c)
+	}
 
 	// Set Q = I.
 	for i := 0; i < c; i++ {
@@ -132,6 +144,8 @@ func (lq *LQ) QTo(dst *Dense) {
 			1, h, qCopy.mat,
 			0, dst.mat)
 	}
+
+	return dst
 }
 
 // SolveLQ finds a minimum-norm solution to a system of linear equations defined
diff --git a/matrix/mat64/lq_test.go b/matrix/mat64/lq_test.go
index 657b8b48..746283b4 100644
--- a/matrix/mat64/lq_test.go
+++ b/matrix/mat64/lq_test.go
@@ -29,17 +29,16 @@ func TestLQ(t *testing.T) {
 
 		var lq LQ
 		lq.Factorize(a)
-		var l, q Dense
-		lq.QTo(&q)
+		q := lq.QTo(nil)
 
-		if !isOrthonormal(&q, 1e-10) {
+		if !isOrthonormal(q, 1e-10) {
 			t.Errorf("Q is not orthonormal: m = %v, n = %v", m, n)
 		}
 
-		lq.LTo(&l)
+		l := lq.LTo(nil)
 
 		var got Dense
-		got.Mul(&l, &q)
+		got.Mul(l, q)
 		if !EqualApprox(&got, &want, 1e-12) {
 			t.Errorf("LQ does not equal original matrix. \nWant: %v\nGot: %v", want, got)
 		}
diff --git a/matrix/mat64/lu.go b/matrix/mat64/lu.go
index a059a5a7..ffacc742 100644
--- a/matrix/mat64/lu.go
+++ b/matrix/mat64/lu.go
@@ -205,9 +205,14 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y *Vector) {
 }
 
 // LTo extracts the lower triangular matrix from an LU factorization.
-func (lu *LU) LTo(dst *TriDense) {
+// If dst is nil, a new matrix is allocated. The resulting L matrix is returned.
+func (lu *LU) LTo(dst *TriDense) *TriDense {
 	_, n := lu.lu.Dims()
-	dst.reuseAs(n, false)
+	if dst == nil {
+		dst = NewTriDense(n, matrix.Lower, nil)
+	} else {
+		dst.reuseAs(n, matrix.Lower)
+	}
 	// Extract the lower triangular elements.
 	for i := 0; i < n; i++ {
 		for j := 0; j < i; j++ {
@@ -218,18 +223,25 @@ func (lu *LU) LTo(dst *TriDense) {
 	for i := 0; i < n; i++ {
 		dst.mat.Data[i*dst.mat.Stride+i] = 1
 	}
+	return dst
 }
 
 // UTo extracts the upper triangular matrix from an LU factorization.
-func (lu *LU) UTo(dst *TriDense) {
+// If dst is nil, a new matrix is allocated. The resulting U matrix is returned.
+func (lu *LU) UTo(dst *TriDense) *TriDense {
 	_, n := lu.lu.Dims()
-	dst.reuseAs(n, true)
+	if dst == nil {
+		dst = NewTriDense(n, matrix.Upper, nil)
+	} else {
+		dst.reuseAs(n, matrix.Upper)
+	}
 	// Extract the upper triangular elements.
 	for i := 0; i < n; i++ {
 		for j := i; j < n; j++ {
 			dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j]
 		}
 	}
+	return dst
 }
 
 // Permutation constructs an r×r permutation matrix with the given row swaps.
diff --git a/matrix/mat64/lu_test.go b/matrix/mat64/lu_test.go
index 6468715e..a2693d21 100644
--- a/matrix/mat64/lu_test.go
+++ b/matrix/mat64/lu_test.go
@@ -23,15 +23,13 @@ func TestLUD(t *testing.T) {
 		var lu LU
 		lu.Factorize(a)
 
-		var l, u TriDense
-		lu.LTo(&l)
-		lu.UTo(&u)
+		l := lu.LTo(nil)
+		u := lu.UTo(nil)
 		var p Dense
 		pivot := lu.Pivot(nil)
 		p.Permutation(n, pivot)
 		var got Dense
-		got.Mul(&p, &l)
-		got.Mul(&got, &u)
+		got.Product(&p, l, u)
 		if !EqualApprox(&got, &want, 1e-12) {
 			t.Errorf("PLU does not equal original matrix.\nWant: %v\n Got: %v", want, got)
 		}
diff --git a/matrix/mat64/qr.go b/matrix/mat64/qr.go
index 770bd3b3..b25bb2af 100644
--- a/matrix/mat64/qr.go
+++ b/matrix/mat64/qr.go
@@ -61,9 +61,14 @@ func (qr *QR) Factorize(a Matrix) {
 // and upper triangular matrices.
 
 // RTo extracts the m×n upper trapezoidal matrix from a QR decomposition.
-func (qr *QR) RTo(dst *Dense) {
+// If dst is nil, a new matrix is allocated. The resulting dst matrix is returned.
+func (qr *QR) RTo(dst *Dense) *Dense {
 	r, c := qr.qr.Dims()
-	dst.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	// Disguise the QR as an upper triangular
 	t := &TriDense{
@@ -82,12 +87,19 @@ func (qr *QR) RTo(dst *Dense) {
 	for i := r; i < c; i++ {
 		zero(dst.mat.Data[i*dst.mat.Stride : i*dst.mat.Stride+c])
 	}
+
+	return dst
 }
 
 // QTo extracts the m×m orthonormal matrix Q from a QR decomposition.
-func (qr *QR) QTo(dst *Dense) {
+// If dst is nil, a new matrix is allocated. The resulting Q matrix is returned.
+func (qr *QR) QTo(dst *Dense) *Dense {
 	r, _ := qr.qr.Dims()
-	dst.reuseAsZeroed(r, r)
+	if dst == nil {
+		dst = NewDense(r, r, nil)
+	} else {
+		dst.reuseAsZeroed(r, r)
+	}
 
 	// Set Q = I.
 	for i := 0; i < r*r; i += r + 1 {
@@ -99,6 +111,8 @@ func (qr *QR) QTo(dst *Dense) {
 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, -1)
 	work = make([]float64, int(work[0]))
 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, len(work))
+
+	return dst
 }
 
 // SolveQR finds a minimum-norm solution to a system of linear equations defined
diff --git a/matrix/mat64/qr_test.go b/matrix/mat64/qr_test.go
index 61407137..6914c7e8 100644
--- a/matrix/mat64/qr_test.go
+++ b/matrix/mat64/qr_test.go
@@ -32,17 +32,16 @@ func TestQR(t *testing.T) {
 
 		var qr QR
 		qr.Factorize(a)
-		var q, r Dense
-		qr.QTo(&q)
+		q := qr.QTo(nil)
 
-		if !isOrthonormal(&q, 1e-10) {
+		if !isOrthonormal(q, 1e-10) {
 			t.Errorf("Q is not orthonormal: m = %v, n = %v", m, n)
 		}
 
-		qr.RTo(&r)
+		r := qr.RTo(nil)
 
 		var got Dense
-		got.Mul(&q, &r)
+		got.Mul(q, r)
 		if !EqualApprox(&got, &want, 1e-12) {
 			t.Errorf("QR does not equal original matrix. \nWant: %v\nGot: %v", want, got)
 		}
diff --git a/matrix/mat64/svd.go b/matrix/mat64/svd.go
index bd8f4992..784e5d9d 100644
--- a/matrix/mat64/svd.go
+++ b/matrix/mat64/svd.go
@@ -141,14 +141,18 @@ func (svd *SVD) Values(s []float64) []float64 {
 // UTo extracts the matrix U from the singular value decomposition, storing
 // the result in-place into dst. U is size m×m if svd.Kind() == SVDFull,
 // of size m×min(m,n) if svd.Kind() == SVDThin, and UTo panics otherwise.
-func (svd *SVD) UTo(dst *Dense) {
+func (svd *SVD) UTo(dst *Dense) *Dense {
 	kind := svd.kind
 	if kind != matrix.SVDFull && kind != matrix.SVDThin {
 		panic("mat64: improper SVD kind")
 	}
 	r := svd.u.Rows
 	c := svd.u.Cols
-	dst.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     svd.u,
@@ -156,19 +160,25 @@ func (svd *SVD) UTo(dst *Dense) {
 		capCols: c,
 	}
 	dst.Copy(tmp)
+
+	return dst
 }
 
 // VTo extracts the matrix V from the singular value decomposition, storing
 // the result in-place into dst. V is size n×n if svd.Kind() == SVDFull,
 // of size n×min(m,n) if svd.Kind() == SVDThin, and VTo panics otherwise.
-func (svd *SVD) VTo(dst *Dense) {
+func (svd *SVD) VTo(dst *Dense) *Dense {
 	kind := svd.kind
 	if kind != matrix.SVDFull && kind != matrix.SVDThin {
 		panic("mat64: improper SVD kind")
 	}
 	r := svd.vt.Rows
 	c := svd.vt.Cols
-	dst.reuseAs(c, r)
+	if dst == nil {
+		dst = NewDense(c, r, nil)
+	} else {
+		dst.reuseAs(c, r)
+	}
 
 	tmp := &Dense{
 		mat:     svd.vt,
@@ -176,4 +186,6 @@ func (svd *SVD) VTo(dst *Dense) {
 		capCols: c,
 	}
 	dst.Copy(tmp.T())
+
+	return dst
 }
diff --git a/matrix/mat64/svd_test.go b/matrix/mat64/svd_test.go
index 7f1e1b96..34e21a6f 100644
--- a/matrix/mat64/svd_test.go
+++ b/matrix/mat64/svd_test.go
@@ -169,9 +169,5 @@ func TestSVD(t *testing.T) {
 }
 
 func extractSVD(svd *SVD) (s []float64, u, v *Dense) {
-	var um, vm Dense
-	svd.UTo(&um)
-	svd.VTo(&vm)
-	s = svd.Values(nil)
-	return s, &um, &vm
+	return svd.Values(nil), svd.UTo(nil), svd.VTo(nil)
 }