diff --git a/matrix/mat64/cholesky.go b/matrix/mat64/cholesky.go
index 4175a51b..dca84072 100644
--- a/matrix/mat64/cholesky.go
+++ b/matrix/mat64/cholesky.go
@@ -39,7 +39,8 @@ type Cholesky struct {
 // the norm is estimated from the decomposition.
 func (c *Cholesky) updateCond(norm float64) {
 	n := c.chol.mat.N
-	work := make([]float64, 3*n)
+	work := getFloats(3*n, false)
+	defer putFloats(work)
 	if norm < 0 {
 		// This is an approximation. By the definition of a norm, ||AB|| <= ||A|| ||B||.
 		// Here, A = U^T * U.
@@ -52,8 +53,9 @@ func (c *Cholesky) updateCond(norm float64) {
 		norm = unorm * lnorm
 	}
 	sym := c.chol.asSymBlas()
-	iwork := make([]int, n)
+	iwork := getInts(n, false)
 	v := lapack64.Pocon(sym, norm, work, iwork)
+	putInts(iwork)
 	c.cond = 1 / v
 }
 
@@ -70,8 +72,9 @@ func (c *Cholesky) Factorize(a Symmetric) (ok bool) {
 	copySymIntoTriangle(c.chol, a)
 
 	sym := c.chol.asSymBlas()
-	work := make([]float64, c.chol.mat.N)
+	work := getFloats(c.chol.mat.N, false)
 	norm := lapack64.Lansy(matrix.CondNorm, sym, work)
+	putFloats(work)
 	_, ok = lapack64.Potrf(sym)
 	if ok {
 		c.updateCond(norm)
@@ -345,7 +348,8 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x *Vector) (ok bool
 	//   EPFL Technical Report 161468 (2004)
 	//   http://infoscience.epfl.ch/record/161468
 
-	work := make([]float64, n)
+	work := getFloats(n, false)
+	defer putFloats(work)
 	blas64.Copy(n, x.RawVector(), blas64.Vector{1, work})
 
 	if alpha > 0 {
@@ -404,8 +408,10 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x *Vector) (ok bool
 		return false
 	}
 	norm = math.Sqrt((1 + norm) * (1 - norm))
-	cos := make([]float64, n)
-	sin := make([]float64, n)
+	cos := getFloats(n, false)
+	defer putFloats(cos)
+	sin := getFloats(n, false)
+	defer putFloats(sin)
 	for i := n - 1; i >= 0; i-- {
 		// Compute parameters of Givens matrices that zero elements of p
 		// backwards.
diff --git a/matrix/mat64/dense_arithmetic.go b/matrix/mat64/dense_arithmetic.go
index 8b94bbb9..bc7945f3 100644
--- a/matrix/mat64/dense_arithmetic.go
+++ b/matrix/mat64/dense_arithmetic.go
@@ -230,18 +230,22 @@ func (m *Dense) Inverse(a Matrix) error {
 	default:
 		m.Copy(a)
 	}
-	ipiv := make([]int, r)
+	ipiv := getInts(r, false)
+	defer putInts(ipiv)
 	ok := lapack64.Getrf(m.mat, ipiv)
 	if !ok {
 		return matrix.Condition(math.Inf(1))
 	}
-	work := make([]float64, 1, 4*r) // must be at least 4*r for cond.
+	work := getFloats(4*r, false) // must be at least 4*r for cond.
 	lapack64.Getri(m.mat, ipiv, work, -1)
 	if int(work[0]) > 4*r {
-		work = make([]float64, int(work[0]))
+		l := int(work[0])
+		putFloats(work)
+		work = getFloats(l, false)
 	} else {
 		work = work[:4*r]
 	}
+	defer putFloats(work)
 	lapack64.Getri(m.mat, ipiv, work, len(work))
 	norm := lapack64.Lange(matrix.CondNorm, m.mat, work)
 	rcond := lapack64.Gecon(matrix.CondNorm, m.mat, norm, work, ipiv) // reuse ipiv
@@ -426,7 +430,8 @@ func (m *Dense) Mul(a, b Matrix) {
 		}
 	}
 
-	row := make([]float64, ac)
+	row := getFloats(ac, false)
+	defer putFloats(row)
 	for r := 0; r < ar; r++ {
 		for i := range row {
 			row[i] = a.At(r, i)
diff --git a/matrix/mat64/eigen.go b/matrix/mat64/eigen.go
index 81fc9fc4..08652864 100644
--- a/matrix/mat64/eigen.go
+++ b/matrix/mat64/eigen.go
@@ -44,11 +44,12 @@ func (e *EigenSym) Factorize(a Symmetric, vectors bool) (ok bool) {
 		jobz = lapack.ComputeEV
 	}
 	w := make([]float64, n)
-	work := make([]float64, 1)
+	work := []float64{0}
 	lapack64.Syev(jobz, sd.mat, w, work, -1)
 
-	work = make([]float64, int(work[0]))
+	work = getFloats(int(work[0]), false)
 	ok = lapack64.Syev(jobz, sd.mat, w, work, len(work))
+	putFloats(work)
 	if !ok {
 		e.vectorsComputed = false
 		e.values = nil
@@ -165,13 +166,16 @@ func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) {
 		jobvr = lapack.ComputeRightEV
 	}
 
-	wr := make([]float64, c)
-	wi := make([]float64, c)
+	wr := getFloats(c, false)
+	defer putFloats(wr)
+	wi := getFloats(c, false)
+	defer putFloats(wi)
 
-	work := make([]float64, 1)
+	work := []float64{0}
 	lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, -1)
-	work = make([]float64, int(work[0]))
+	work = getFloats(int(work[0]), false)
 	first := lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, len(work))
+	putFloats(work)
 
 	if first != 0 {
 		e.values = nil
diff --git a/matrix/mat64/lq.go b/matrix/mat64/lq.go
index 61892d6d..30cf5d02 100644
--- a/matrix/mat64/lq.go
+++ b/matrix/mat64/lq.go
@@ -24,11 +24,13 @@ func (lq *LQ) updateCond() {
 	// A = LQ, where Q is orthonormal. Orthonormal multiplications do not change
 	// the condition number. Thus, ||A|| = ||L|| ||Q|| = ||Q||.
 	m := lq.lq.mat.Rows
-	work := make([]float64, 3*m)
-	iwork := make([]int, m)
+	work := getFloats(3*m, false)
+	iwork := getInts(m, false)
 	l := lq.lq.asTriDense(m, blas.NonUnit, blas.Lower)
 	v := lapack64.Trcon(matrix.CondNorm, l.mat, work, iwork)
 	lq.cond = 1 / v
+	putFloats(work)
+	putInts(iwork)
 }
 
 // Factorize computes the LQ factorization of an m×n matrix a where n <= m. The LQ
@@ -47,11 +49,12 @@ func (lq *LQ) Factorize(a Matrix) {
 		lq.lq = &Dense{}
 	}
 	lq.lq.Clone(a)
-	work := make([]float64, 1)
+	work := []float64{0}
 	lq.tau = make([]float64, k)
 	lapack64.Gelqf(lq.lq.mat, lq.tau, work, -1)
-	work = make([]float64, int(work[0]))
+	work = getFloats(int(work[0]), false)
 	lapack64.Gelqf(lq.lq.mat, lq.tau, work, len(work))
+	putFloats(work)
 	lq.updateCond()
 }
 
@@ -110,10 +113,11 @@ func (lq *LQ) QTo(dst *Dense) *Dense {
 	}
 
 	// Construct Q from the elementary reflectors.
-	work := make([]float64, 1)
+	work := []float64{0}
 	lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, -1)
-	work = make([]float64, int(work[0]))
+	work = getFloats(int(work[0]), false)
 	lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, len(work))
+	putFloats(work)
 
 	return dst
 }
@@ -152,10 +156,11 @@ func (m *Dense) SolveLQ(lq *LQ, trans bool, b Matrix) error {
 	x.Copy(b)
 	t := lq.lq.asTriDense(lq.lq.mat.Rows, blas.NonUnit, blas.Lower).mat
 	if trans {
-		work := make([]float64, 1)
+		work := []float64{0}
 		lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, x.mat, work, -1)
-		work = make([]float64, int(work[0]))
+		work = getFloats(int(work[0]), false)
 		lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, x.mat, work, len(work))
+		putFloats(work)
 
 		ok := lapack64.Trtrs(blas.Trans, t, x.mat)
 		if !ok {
@@ -169,10 +174,11 @@ func (m *Dense) SolveLQ(lq *LQ, trans bool, b Matrix) error {
 		for i := r; i < c; i++ {
 			zero(x.mat.Data[i*x.mat.Stride : i*x.mat.Stride+bc])
 		}
-		work := make([]float64, 1)
+		work := []float64{0}
 		lapack64.Ormlq(blas.Left, blas.Trans, lq.lq.mat, lq.tau, x.mat, work, -1)
-		work = make([]float64, int(work[0]))
+		work = getFloats(int(work[0]), false)
 		lapack64.Ormlq(blas.Left, blas.Trans, lq.lq.mat, lq.tau, x.mat, work, len(work))
+		putFloats(work)
 	}
 	// M was set above to be the correct size for the result.
 	m.Copy(x)
diff --git a/matrix/mat64/lu.go b/matrix/mat64/lu.go
index ffacc742..6db840f8 100644
--- a/matrix/mat64/lu.go
+++ b/matrix/mat64/lu.go
@@ -27,8 +27,10 @@ type LU struct {
 // norm of the original matrix. If norm is negative it will be estimated.
 func (lu *LU) updateCond(norm float64) {
 	n := lu.lu.mat.Cols
-	work := make([]float64, 4*n)
-	iwork := make([]int, n)
+	work := getFloats(4*n, false)
+	defer putFloats(work)
+	iwork := getInts(n, false)
+	defer putInts(iwork)
 	if norm < 0 {
 		// This is an approximation. By the definition of a norm, ||AB|| <= ||A|| ||B||.
 		// The condition number is ||A|| || A^-1||, so this will underestimate
@@ -68,8 +70,9 @@ func (lu *LU) Factorize(a Matrix) {
 		lu.pivot = make([]int, r)
 	}
 	lu.pivot = lu.pivot[:r]
-	work := make([]float64, r)
+	work := getFloats(r, false)
 	anorm := lapack64.Lange(matrix.CondNorm, lu.lu.mat, work)
+	putFloats(work)
 	lapack64.Getrf(lu.lu.mat, lu.pivot)
 	lu.updateCond(anorm)
 }
@@ -99,7 +102,8 @@ func (lu *LU) Det() float64 {
 // division expressions is generally improved by working in log space.
 func (lu *LU) LogDet() (det float64, sign float64) {
 	_, n := lu.lu.Dims()
-	logDiag := make([]float64, n)
+	logDiag := getFloats(n, false)
+	defer putFloats(logDiag)
 	sign = 1.0
 	for i := 0; i < n; i++ {
 		v := lu.lu.at(i, i)
@@ -171,8 +175,10 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y *Vector) {
 		lu.lu.Copy(orig.lu)
 	}
 
-	xs := make([]float64, n)
-	ys := make([]float64, n)
+	xs := getFloats(n, false)
+	defer putFloats(xs)
+	ys := getFloats(n, false)
+	defer putFloats(ys)
 	for i := 0; i < n; i++ {
 		xs[i] = x.at(i)
 		ys[i] = y.at(i)
diff --git a/matrix/mat64/matrix.go b/matrix/mat64/matrix.go
index 54d87d92..0b5d3714 100644
--- a/matrix/mat64/matrix.go
+++ b/matrix/mat64/matrix.go
@@ -292,48 +292,58 @@ func Cond(a Matrix, norm float64) float64 {
 		// Use the LU decomposition to compute the condition number.
 		tmp := getWorkspace(m, n, false)
 		tmp.Copy(a)
-		work := make([]float64, 4*n)
+		work := getFloats(4*n, false)
 		aNorm := lapack64.Lange(lnorm, tmp.mat, work)
-		pivot := make([]int, m)
+		pivot := getInts(m, false)
 		lapack64.Getrf(tmp.mat, pivot)
 		iwork := make([]int, n)
 		v := lapack64.Gecon(lnorm, tmp.mat, aNorm, work, iwork)
 		putWorkspace(tmp)
+		putFloats(work)
+		putInts(pivot)
 		return 1 / v
 	}
 	if m > n {
 		// Use the QR factorization to compute the condition number.
 		tmp := getWorkspace(m, n, false)
 		tmp.Copy(a)
-		work := make([]float64, 3*n)
-		tau := make([]float64, min(m, n))
+		work := getFloats(3*n, false)
+		tau := getFloats(min(m, n), false)
 		lapack64.Geqrf(tmp.mat, tau, work, -1)
-		if int(work[0]) > len(work) {
-			work = make([]float64, int(work[0]))
+		if l := int(work[0]); l > len(work) {
+			putFloats(work)
+			work = getFloats(l, false)
 		}
 		lapack64.Geqrf(tmp.mat, tau, work, len(work))
 
-		iwork := make([]int, n)
+		iwork := getInts(n, false)
 		r := tmp.asTriDense(n, blas.NonUnit, blas.Upper)
 		v := lapack64.Trcon(lnorm, r.mat, work, iwork)
 		putWorkspace(tmp)
+		putFloats(work)
+		putFloats(tau)
+		putInts(iwork)
 		return 1 / v
 	}
 	// Use the LQ factorization to compute the condition number.
 	tmp := getWorkspace(m, n, false)
 	tmp.Copy(a)
-	work := make([]float64, 3*m)
-	tau := make([]float64, min(m, n))
+	work := getFloats(3*m, false)
+	tau := getFloats(min(m, n), false)
 	lapack64.Gelqf(tmp.mat, tau, work, -1)
-	if int(work[0]) > len(work) {
-		work = make([]float64, int(work[0]))
+	if l := int(work[0]); l > len(work) {
+		putFloats(work)
+		work = getFloats(l, false)
 	}
 	lapack64.Gelqf(tmp.mat, tau, work, len(work))
 
-	iwork := make([]int, m)
+	iwork := getInts(m, false)
 	l := tmp.asTriDense(m, blas.NonUnit, blas.Lower)
 	v := lapack64.Trcon(lnorm, l.mat, work, iwork)
 	putWorkspace(tmp)
+	putFloats(work)
+	putFloats(tau)
+	putInts(iwork)
 	return 1 / v
 }
 
@@ -679,21 +689,24 @@ func Norm(a Matrix, norm float64) float64 {
 		rm := rma.RawMatrix()
 		n := normLapack(norm, aTrans)
 		if n == lapack.MaxColumnSum {
-			work = make([]float64, rm.Cols)
+			work = getFloats(rm.Cols, false)
+			defer putFloats(work)
 		}
 		return lapack64.Lange(n, rm, work)
 	case RawTriangular:
 		rm := rma.RawTriangular()
 		n := normLapack(norm, aTrans)
 		if n == lapack.MaxRowSum || n == lapack.MaxColumnSum {
-			work = make([]float64, rm.N)
+			work = getFloats(rm.N, false)
+			defer putFloats(work)
 		}
 		return lapack64.Lantr(n, rm, work)
 	case RawSymmetricer:
 		rm := rma.RawSymmetric()
 		n := normLapack(norm, aTrans)
 		if n == lapack.MaxRowSum || n == lapack.MaxColumnSum {
-			work = make([]float64, rm.N)
+			work = getFloats(rm.N, false)
+			defer putFloats(work)
 		}
 		return lapack64.Lansy(n, rm, work)
 	case *Vector:
diff --git a/matrix/mat64/pool.go b/matrix/mat64/pool.go
index 53f04df4..10eeef06 100644
--- a/matrix/mat64/pool.go
+++ b/matrix/mat64/pool.go
@@ -54,6 +54,12 @@ var (
 
 	// poolVec is the Vector equivalent of pool.
 	poolVec [63]sync.Pool
+
+	// poolFloats is the []float64 equivalent of pool.
+	poolFloats [63]sync.Pool
+
+	// poolInts is the []int equivalent of pool.
+	poolInts [63]sync.Pool
 )
 
 func init() {
@@ -81,6 +87,12 @@ func init() {
 				Data: make([]float64, l),
 			}}
 		}
+		poolFloats[i].New = func() interface{} {
+			return make([]float64, l)
+		}
+		poolInts[i].New = func() interface{} {
+			return make([]int, l)
+		}
 	}
 }
 
@@ -182,3 +194,41 @@ func getWorkspaceVec(n int, clear bool) *Vector {
 func putWorkspaceVec(v *Vector) {
 	poolVec[bits(uint64(cap(v.mat.Data)))].Put(v)
 }
+
+// getFloats returns a []float64 of length l and a cap that is
+// less than 2*l. If clear is true, the slice visible is zeroed.
+func getFloats(l int, clear bool) []float64 {
+	w := poolFloats[bits(uint64(l))].Get().([]float64)
+	w = w[:l]
+	if clear {
+		zero(w)
+	}
+	return w
+}
+
+// putFloats replaces a used []float64 into the appropriate size
+// workspace pool. putFloats must not be called with a slice
+// where references to the underlying data have been kept.
+func putFloats(w []float64) {
+	poolFloats[bits(uint64(cap(w)))].Put(w)
+}
+
+// getInts returns a []ints of length l and a cap that is
+// less than 2*l. If clear is true, the slice visible is zeroed.
+func getInts(l int, clear bool) []int {
+	w := poolInts[bits(uint64(l))].Get().([]int)
+	w = w[:l]
+	if clear {
+		for i := range w {
+			w[i] = 0
+		}
+	}
+	return w
+}
+
+// putInts replaces a used []int into the appropriate size
+// workspace pool. putInts must not be called with a slice
+// where references to the underlying data have been kept.
+func putInts(w []int) {
+	poolInts[bits(uint64(cap(w)))].Put(w)
+}
diff --git a/matrix/mat64/qr.go b/matrix/mat64/qr.go
index b25bb2af..0dfdb6e3 100644
--- a/matrix/mat64/qr.go
+++ b/matrix/mat64/qr.go
@@ -25,10 +25,12 @@ func (qr *QR) updateCond() {
 	// A = QR, where Q is orthonormal. Orthonormal multiplications do not change
 	// the condition number. Thus, ||A|| = ||Q|| ||R|| = ||R||.
 	n := qr.qr.mat.Cols
-	work := make([]float64, 3*n)
-	iwork := make([]int, n)
+	work := getFloats(3*n, false)
+	iwork := getInts(n, false)
 	r := qr.qr.asTriDense(n, blas.NonUnit, blas.Upper)
 	v := lapack64.Trcon(matrix.CondNorm, r.mat, work, iwork)
+	putFloats(work)
+	putInts(iwork)
 	qr.cond = 1 / v
 }
 
@@ -48,12 +50,13 @@ func (qr *QR) Factorize(a Matrix) {
 		qr.qr = &Dense{}
 	}
 	qr.qr.Clone(a)
-	work := make([]float64, 1)
+	work := []float64{0}
 	qr.tau = make([]float64, k)
 	lapack64.Geqrf(qr.qr.mat, qr.tau, work, -1)
 
-	work = make([]float64, int(work[0]))
+	work = getFloats(int(work[0]), false)
 	lapack64.Geqrf(qr.qr.mat, qr.tau, work, len(work))
+	putFloats(work)
 	qr.updateCond()
 }
 
@@ -107,10 +110,11 @@ func (qr *QR) QTo(dst *Dense) *Dense {
 	}
 
 	// Construct Q from the elementary reflectors.
-	work := make([]float64, 1)
+	work := []float64{0}
 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, -1)
-	work = make([]float64, int(work[0]))
+	work = getFloats(int(work[0]), false)
 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, len(work))
+	putFloats(work)
 
 	return dst
 }
@@ -156,15 +160,17 @@ func (m *Dense) SolveQR(qr *QR, trans bool, b Matrix) error {
 		for i := c; i < r; i++ {
 			zero(x.mat.Data[i*x.mat.Stride : i*x.mat.Stride+bc])
 		}
-		work := make([]float64, 1)
+		work := []float64{0}
 		lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, x.mat, work, -1)
-		work = make([]float64, int(work[0]))
+		work = getFloats(int(work[0]), false)
 		lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, x.mat, work, len(work))
+		putFloats(work)
 	} else {
-		work := make([]float64, 1)
+		work := []float64{0}
 		lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, x.mat, work, -1)
-		work = make([]float64, int(work[0]))
+		work = getFloats(int(work[0]), false)
 		lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, x.mat, work, len(work))
+		putFloats(work)
 
 		ok := lapack64.Trtrs(blas.NoTrans, t, x.mat)
 		if !ok {
diff --git a/matrix/mat64/solve.go b/matrix/mat64/solve.go
index 001ed330..e30bada7 100644
--- a/matrix/mat64/solve.go
+++ b/matrix/mat64/solve.go
@@ -64,9 +64,11 @@ func (m *Dense) Solve(a, b Matrix) error {
 
 		rm := rma.RawTriangular()
 		blas64.Trsm(side, tA, 1, rm, m.mat)
-		work := make([]float64, 3*rm.N)
-		iwork := make([]int, rm.N)
+		work := getFloats(3*rm.N, false)
+		iwork := getInts(rm.N, false)
 		cond := lapack64.Trcon(matrix.CondNorm, rm, work, iwork)
+		putFloats(work)
+		putInts(iwork)
 		if cond > matrix.ConditionTolerance {
 			return matrix.Condition(cond)
 		}
diff --git a/matrix/mat64/svd.go b/matrix/mat64/svd.go
index 784e5d9d..0a4f2062 100644
--- a/matrix/mat64/svd.go
+++ b/matrix/mat64/svd.go
@@ -92,10 +92,11 @@ func (svd *SVD) Factorize(a Matrix, kind matrix.SVDKind) (ok bool) {
 	svd.kind = kind
 	svd.s = use(svd.s, min(m, n))
 
-	work := make([]float64, 1)
+	work := []float64{0}
 	lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, -1)
-	work = make([]float64, int(work[0]))
+	work = getFloats(int(work[0]), false)
 	ok = lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, len(work))
+	putFloats(work)
 	if !ok {
 		svd.kind = 0
 	}
diff --git a/matrix/mat64/triangular.go b/matrix/mat64/triangular.go
index 99f0f6cf..394d5159 100644
--- a/matrix/mat64/triangular.go
+++ b/matrix/mat64/triangular.go
@@ -331,9 +331,11 @@ func (t *TriDense) InverseTri(a Triangular) error {
 	n, _ := a.Triangle()
 	t.reuseAs(a.Triangle())
 	t.Copy(a)
-	work := make([]float64, 3*n)
-	iwork := make([]int, n)
+	work := getFloats(3*n, false)
+	iwork := getInts(n, false)
 	cond := lapack64.Trcon(matrix.CondNorm, t.mat, work, iwork)
+	putFloats(work)
+	putInts(iwork)
 	if math.IsInf(cond, 1) {
 		return matrix.Condition(cond)
 	}