matrix/mat64: use pools for []float64 and []int where possible
MODIFY
matrix/mat64/cholesky.go
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.

MODIFY
matrix/mat64/dense_arithmetic.go
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)

MODIFY
matrix/mat64/eigen.go
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

MODIFY
matrix/mat64/lu.go
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)

MODIFY
matrix/mat64/matrix.go
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:

MODIFY
matrix/mat64/pool.go
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)
+}

MODIFY
matrix/mat64/qr.go
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 {

MODIFY
matrix/mat64/solve.go
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)
}

MODIFY
matrix/mat64/svd.go
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
}

MODIFY
matrix/mat64/triangular.go
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)
}

matrix/mat64: use pools for []float64 and []int where possible