diff --git a/mat64/lq.go b/mat64/lq.go
index b351136..7d2331c 100644
--- a/mat64/lq.go
+++ b/mat64/lq.go
@@ -36,7 +36,7 @@ func (lq *LQ) updateCond() {
 //
 // The LQ decomposition is a factorization of the matrix A such that A = L * Q.
 // The matrix Q is an orthonormal n×n matrix, and L is an m×n upper triangular matrix.
-// L and Q can be extracted from the LFromLQ and QFromLQ methods on Dense.
+// L and Q can be extracted from the LTo and QTo methods.
 func (lq *LQ) Factorize(a Matrix) {
 	m, n := a.Dims()
 	if m > n {
@@ -58,10 +58,15 @@ func (lq *LQ) Factorize(a Matrix) {
 // TODO(btracey): Add in the "Reduced" forms for extracting the m×m orthogonal
 // and upper triangular matrices.
 
-// LFromLQ extracts the m×n lower trapezoidal matrix from a LQ decomposition.
-func (m *Dense) LFromLQ(lq *LQ) {
+// LTo extracts the m×n lower trapezoidal matrix from a LQ decomposition.
+// 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()
-	m.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	// Disguise the LQ as a lower triangular
 	t := &TriDense{
@@ -74,25 +79,32 @@ func (m *Dense) LFromLQ(lq *LQ) {
 		},
 		cap: lq.lq.capCols,
 	}
-	m.Copy(t)
+	dst.Copy(t)
 
 	if r == c {
-		return
+		return dst
 	}
 	// Zero right of the triangular.
 	for i := 0; i < r; i++ {
-		zero(m.mat.Data[i*m.mat.Stride+r : i*m.mat.Stride+c])
+		zero(dst.mat.Data[i*dst.mat.Stride+r : i*dst.mat.Stride+c])
 	}
+
+	return dst
 }
 
-// QFromLQ extracts the n×n orthonormal matrix Q from an LQ decomposition.
-func (m *Dense) QFromLQ(lq *LQ) {
+// QTo extracts the n×n orthonormal matrix Q from an LQ decomposition.
+// 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()
-	m.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++ {
-		v := m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c]
+		v := dst.mat.Data[i*dst.mat.Stride : i*dst.mat.Stride+c]
 		zero(v)
 		v[i] = 1
 	}
@@ -127,11 +139,13 @@ func (m *Dense) QFromLQ(lq *LQ) {
 
 		// Compute the multiplication matrix.
 		blas64.Ger(-lq.tau[i], v, v, h)
-		qCopy.Copy(m)
+		qCopy.Copy(dst)
 		blas64.Gemm(blas.NoTrans, blas.NoTrans,
 			1, h, qCopy.mat,
-			0, m.mat)
+			0, dst.mat)
 	}
+
+	return dst
 }
 
 // SolveLQ finds a minimum-norm solution to a system of linear equations defined
diff --git a/mat64/lq_test.go b/mat64/lq_test.go
index 05478ce..746283b 100644
--- a/mat64/lq_test.go
+++ b/mat64/lq_test.go
@@ -27,19 +27,18 @@ func TestLQ(t *testing.T) {
 		var want Dense
 		want.Clone(a)
 
-		lq := &LQ{}
+		var lq LQ
 		lq.Factorize(a)
-		var l, q Dense
-		q.QFromLQ(lq)
+		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)
 		}
 
-		l.LFromLQ(lq)
+		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/mat64/lu.go b/mat64/lu.go
index 6115b00..603f259 100644
--- a/mat64/lu.go
+++ b/mat64/lu.go
@@ -51,7 +51,7 @@ func (lu *LU) updateCond(norm float64) {
 // The LU factorization is computed with pivoting, and so really the decomposition
 // is a PLU decomposition where P is a permutation matrix. The individual matrix
 // factors can be extracted from the factorization using the Permutation method
-// on Dense, and the LFrom and UFrom methods on TriDense.
+// on Dense, and the LU LTo and UTo methods.
 func (lu *LU) Factorize(a Matrix) {
 	r, c := a.Dims()
 	if r != c {
@@ -204,32 +204,44 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y *Vector) {
 	lu.updateCond(-1)
 }
 
-// LFromLU extracts the lower triangular matrix from an LU factorization.
-func (t *TriDense) LFromLU(lu *LU) {
+// LTo extracts the lower triangular matrix from an LU factorization.
+// 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()
-	t.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++ {
-			t.mat.Data[i*t.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j]
+			dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j]
 		}
 	}
 	// Set ones on the diagonal.
 	for i := 0; i < n; i++ {
-		t.mat.Data[i*t.mat.Stride+i] = 1
+		dst.mat.Data[i*dst.mat.Stride+i] = 1
 	}
+	return dst
 }
 
-// UFromLU extracts the upper triangular matrix from an LU factorization.
-func (t *TriDense) UFromLU(lu *LU) {
+// UTo extracts the upper triangular matrix from an LU factorization.
+// 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()
-	t.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++ {
-			t.mat.Data[i*t.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+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/mat64/lu_test.go b/mat64/lu_test.go
index ed7dd6d..a2693d2 100644
--- a/mat64/lu_test.go
+++ b/mat64/lu_test.go
@@ -20,18 +20,16 @@ func TestLUD(t *testing.T) {
 		var want Dense
 		want.Clone(a)
 
-		lu := &LU{}
+		var lu LU
 		lu.Factorize(a)
 
-		var l, u TriDense
-		l.LFromLU(lu)
-		u.UFromLU(lu)
+		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)
 		}
@@ -93,8 +91,8 @@ func TestLURankOne(t *testing.T) {
 // luReconstruct reconstructs the original A matrix from an LU decomposition.
 func luReconstruct(lu *LU) *Dense {
 	var L, U TriDense
-	L.LFromLU(lu)
-	U.UFromLU(lu)
+	lu.LTo(&L)
+	lu.UTo(&U)
 	var P Dense
 	pivot := lu.Pivot(nil)
 	P.Permutation(len(pivot), pivot)
diff --git a/mat64/qr.go b/mat64/qr.go
index 140f7dc..d9781ec 100644
--- a/mat64/qr.go
+++ b/mat64/qr.go
@@ -37,7 +37,7 @@ func (qr *QR) updateCond() {
 //
 // The QR decomposition is a factorization of the matrix A such that A = Q * R.
 // The matrix Q is an orthonormal m×m matrix, and R is an m×n upper triangular matrix.
-// Q and R can be extracted from the QFromQR and RFromQR methods on Dense.
+// Q and R can be extracted using the QTo and RTo methods.
 func (qr *QR) Factorize(a Matrix) {
 	m, n := a.Dims()
 	if m < n {
@@ -60,10 +60,15 @@ func (qr *QR) Factorize(a Matrix) {
 // TODO(btracey): Add in the "Reduced" forms for extracting the n×n orthogonal
 // and upper triangular matrices.
 
-// RFromQR extracts the m×n upper trapezoidal matrix from a QR decomposition.
-func (m *Dense) RFromQR(qr *QR) {
+// RTo extracts the m×n upper trapezoidal matrix from a QR decomposition.
+// 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()
-	m.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{
@@ -76,29 +81,38 @@ func (m *Dense) RFromQR(qr *QR) {
 		},
 		cap: qr.qr.capCols,
 	}
-	m.Copy(t)
+	dst.Copy(t)
 
 	// Zero below the triangular.
 	for i := r; i < c; i++ {
-		zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c])
+		zero(dst.mat.Data[i*dst.mat.Stride : i*dst.mat.Stride+c])
 	}
+
+	return dst
 }
 
-// QFromQR extracts the m×m orthonormal matrix Q from a QR decomposition.
-func (m *Dense) QFromQR(qr *QR) {
+// QTo extracts the m×m orthonormal matrix Q from a QR decomposition.
+// 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()
-	m.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 {
-		m.mat.Data[i] = 1
+		dst.mat.Data[i] = 1
 	}
 
 	// Construct Q from the elementary reflectors.
 	work := make([]float64, 1)
-	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, m.mat, work, -1)
+	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, m.mat, work, len(work))
+	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/mat64/qr_test.go b/mat64/qr_test.go
index 37a9577..3d58e40 100644
--- a/mat64/qr_test.go
+++ b/mat64/qr_test.go
@@ -30,19 +30,18 @@ func TestQR(t *testing.T) {
 		var want Dense
 		want.Clone(a)
 
-		qr := &QR{}
+		var qr QR
 		qr.Factorize(a)
-		var q, r Dense
-		q.QFromQR(qr)
+		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)
 		}
 
-		r.RFromQR(qr)
+		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/mat64/svd.go b/mat64/svd.go
index 731b48a..f8b99aa 100644
--- a/mat64/svd.go
+++ b/mat64/svd.go
@@ -138,42 +138,54 @@ func (svd *SVD) Values(s []float64) []float64 {
 	return s
 }
 
-// UFromSVD extracts the matrix U from the singular value decomposition, storing
+// UTo extracts the matrix U from the singular value decomposition, storing
 // the result in-place into the receiver. U is size m×m if svd.Kind() == SVDFull,
-// of size m×min(m,n) if svd.Kind() == SVDThin, and UFromSVD panics otherwise.
-func (m *Dense) UFromSVD(svd *SVD) {
+// of size m×min(m,n) if svd.Kind() == SVDThin, and UTo panics otherwise.
+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
-	m.reuseAs(r, c)
+	if dst == nil {
+		dst = NewDense(r, c, nil)
+	} else {
+		dst.reuseAs(r, c)
+	}
 
 	tmp := &Dense{
 		mat:     svd.u,
 		capRows: r,
 		capCols: c,
 	}
-	m.Copy(tmp)
+	dst.Copy(tmp)
+
+	return dst
 }
 
-// VFromSVD extracts the matrix V from the singular value decomposition, storing
+// VTo extracts the matrix V from the singular value decomposition, storing
 // the result in-place into the receiver. V is size n×n if svd.Kind() == SVDFull,
-// of size n×min(m,n) if svd.Kind() == SVDThin, and VFromSVD panics otherwise.
-func (m *Dense) VFromSVD(svd *SVD) {
+// of size n×min(m,n) if svd.Kind() == SVDThin, and VTo panics otherwise.
+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
-	m.reuseAs(c, r)
+	if dst == nil {
+		dst = NewDense(c, r, nil)
+	} else {
+		dst.reuseAs(c, r)
+	}
 
 	tmp := &Dense{
 		mat:     svd.vt,
 		capRows: r,
 		capCols: c,
 	}
-	m.Copy(tmp.T())
+	dst.Copy(tmp.T())
+
+	return dst
 }
diff --git a/mat64/svd_test.go b/mat64/svd_test.go
index 912c110..83e1755 100644
--- a/mat64/svd_test.go
+++ b/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
-	um.UFromSVD(svd)
-	vm.VFromSVD(svd)
-	s = svd.Values(nil)
-	return s, &um, &vm
+	return svd.Values(nil), svd.UTo(nil), svd.VTo(nil)
 }