mat: rename pool helpers to be consistent with type names

This commit is contained in:
Dan Kortschak
2021-05-16 07:40:25 +09:30
parent 608b72c3f0
commit efbee9bf28
21 changed files with 207 additions and 207 deletions

View File

@@ -347,15 +347,15 @@ func (b *BandDense) MulVecTo(dst *VecDense, trans bool, x Vector) {
dst.checkOverlap(xVec.mat) dst.checkOverlap(xVec.mat)
blas64.Gbmv(t, 1, b.mat, xVec.mat, 0, dst.mat) blas64.Gbmv(t, 1, b.mat, xVec.mat, 0, dst.mat)
} else { } else {
xCopy := getWorkspaceVec(n, false) xCopy := getVecDenseWorkspace(n, false)
xCopy.CloneFromVec(xVec) xCopy.CloneFromVec(xVec)
blas64.Gbmv(t, 1, b.mat, xCopy.mat, 0, dst.mat) blas64.Gbmv(t, 1, b.mat, xCopy.mat, 0, dst.mat)
putWorkspaceVec(xCopy) putVecDenseWorkspace(xCopy)
} }
} else { } else {
xCopy := getWorkspaceVec(n, false) xCopy := getVecDenseWorkspace(n, false)
xCopy.CloneFromVec(x) xCopy.CloneFromVec(x)
blas64.Gbmv(t, 1, b.mat, xCopy.mat, 0, dst.mat) blas64.Gbmv(t, 1, b.mat, xCopy.mat, 0, dst.mat)
putWorkspaceVec(xCopy) putVecDenseWorkspace(xCopy)
} }
} }

View File

@@ -226,10 +226,10 @@ func (m *CDense) isolatedWorkspace(a CMatrix) (w *CDense, restore func()) {
if r == 0 || c == 0 { if r == 0 || c == 0 {
panic(ErrZeroLength) panic(ErrZeroLength)
} }
w = getWorkspaceCmplx(r, c, false) w = getCDenseWorkspace(r, c, false)
return w, func() { return w, func() {
m.Copy(w) m.Copy(w)
putWorkspaceCmplx(w) putCDenseWorkspace(w)
} }
} }

View File

@@ -55,8 +55,8 @@ type Cholesky struct {
// the norm is estimated from the decomposition. // the norm is estimated from the decomposition.
func (c *Cholesky) updateCond(norm float64) { func (c *Cholesky) updateCond(norm float64) {
n := c.chol.mat.N n := c.chol.mat.N
work := getFloats(3*n, false) work := getFloat64s(3*n, false)
defer putFloats(work) defer putFloat64s(work)
if norm < 0 { if norm < 0 {
// This is an approximation. By the definition of a norm, // This is an approximation. By the definition of a norm,
// |AB| <= |A| |B|. // |AB| <= |A| |B|.
@@ -139,9 +139,9 @@ func (c *Cholesky) Factorize(a Symmetric) (ok bool) {
copySymIntoTriangle(c.chol, a) copySymIntoTriangle(c.chol, a)
sym := c.chol.asSymBlas() sym := c.chol.asSymBlas()
work := getFloats(c.chol.mat.N, false) work := getFloat64s(c.chol.mat.N, false)
norm := lapack64.Lansy(CondNorm, sym, work) norm := lapack64.Lansy(CondNorm, sym, work)
putFloats(work) putFloat64s(work)
_, ok = lapack64.Potrf(sym) _, ok = lapack64.Potrf(sym)
if ok { if ok {
c.updateCond(norm) c.updateCond(norm)
@@ -592,8 +592,8 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
// EPFL Technical Report 161468 (2004) // EPFL Technical Report 161468 (2004)
// http://infoscience.epfl.ch/record/161468 // http://infoscience.epfl.ch/record/161468
work := getFloats(n, false) work := getFloat64s(n, false)
defer putFloats(work) defer putFloat64s(work)
var xmat blas64.Vector var xmat blas64.Vector
if rv, ok := x.(RawVectorer); ok { if rv, ok := x.(RawVectorer); ok {
xmat = rv.RawVector() xmat = rv.RawVector()
@@ -660,10 +660,10 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
return false return false
} }
norm = math.Sqrt((1 + norm) * (1 - norm)) norm = math.Sqrt((1 + norm) * (1 - norm))
cos := getFloats(n, false) cos := getFloat64s(n, false)
defer putFloats(cos) defer putFloat64s(cos)
sin := getFloats(n, false) sin := getFloat64s(n, false)
defer putFloats(sin) defer putFloat64s(sin)
for i := n - 1; i >= 0; i-- { for i := n - 1; i >= 0; i-- {
// Compute parameters of Givens matrices that zero elements of p // Compute parameters of Givens matrices that zero elements of p
// backwards. // backwards.
@@ -674,8 +674,8 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
sin[i] *= -1 sin[i] *= -1
} }
} }
workMat := getWorkspaceTri(c.chol.mat.N, c.chol.triKind(), false) workMat := getTriDenseWorkspace(c.chol.mat.N, c.chol.triKind(), false)
defer putWorkspaceTri(workMat) defer putTriWorkspace(workMat)
workMat.Copy(c.chol) workMat.Copy(c.chol)
umat := workMat.mat umat := workMat.mat
stride := workMat.mat.Stride stride := workMat.mat.Stride
@@ -752,12 +752,12 @@ func (ch *BandCholesky) Factorize(a SymBanded) (ok bool) {
ch.Reset() ch.Reset()
return false return false
} }
work := getFloats(3*n, false) work := getFloat64s(3*n, false)
iwork := getInts(n, false) iwork := getInts(n, false)
aNorm := lapack64.Lansb(CondNorm, cSym, work) aNorm := lapack64.Lansb(CondNorm, cSym, work)
ch.cond = 1 / lapack64.Pbcon(cSym, aNorm, work, iwork) ch.cond = 1 / lapack64.Pbcon(cSym, aNorm, work, iwork)
putInts(iwork) putInts(iwork)
putFloats(work) putFloat64s(work)
return true return true
} }

View File

@@ -166,10 +166,10 @@ func (m *Dense) isolatedWorkspace(a Matrix) (w *Dense, restore func()) {
if r == 0 || c == 0 { if r == 0 || c == 0 {
panic(ErrZeroLength) panic(ErrZeroLength)
} }
w = getWorkspace(r, c, false) w = getDenseWorkspace(r, c, false)
return w, func() { return w, func() {
m.Copy(w) m.Copy(w)
putWorkspace(w) putDenseWorkspace(w)
} }
} }
@@ -606,8 +606,8 @@ func (m *Dense) Norm(norm float64) float64 {
} }
lnorm := normLapack(norm, false) lnorm := normLapack(norm, false)
if lnorm == lapack.MaxColumnSum { if lnorm == lapack.MaxColumnSum {
work := getFloats(m.mat.Cols, false) work := getFloat64s(m.mat.Cols, false)
defer putFloats(work) defer putFloat64s(work)
return lapack64.Lange(lnorm, m.mat, work) return lapack64.Lange(lnorm, m.mat, work)
} }
return lapack64.Lange(lnorm, m.mat, nil) return lapack64.Lange(lnorm, m.mat, nil)

View File

@@ -226,10 +226,10 @@ func (m *Dense) Inverse(a Matrix) error {
case *Dense: case *Dense:
if m != aU || aTrans { if m != aU || aTrans {
if m == aU || m.checkOverlap(rm.mat) { if m == aU || m.checkOverlap(rm.mat) {
tmp := getWorkspace(r, c, false) tmp := getDenseWorkspace(r, c, false)
tmp.Copy(a) tmp.Copy(a)
m.Copy(tmp) m.Copy(tmp)
putWorkspace(tmp) putDenseWorkspace(tmp)
break break
} }
m.Copy(a) m.Copy(a)
@@ -238,7 +238,7 @@ func (m *Dense) Inverse(a Matrix) error {
m.Copy(a) m.Copy(a)
} }
// Compute the norm of A. // Compute the norm of A.
work := getFloats(4*r, false) // Length must be at least 4*r for Gecon. work := getFloat64s(4*r, false) // Length must be at least 4*r for Gecon.
norm := lapack64.Lange(CondNorm, m.mat, work) norm := lapack64.Lange(CondNorm, m.mat, work)
// Compute the LU factorization of A. // Compute the LU factorization of A.
ipiv := getInts(r, false) ipiv := getInts(r, false)
@@ -256,10 +256,10 @@ func (m *Dense) Inverse(a Matrix) error {
lapack64.Getri(m.mat, ipiv, work, -1) lapack64.Getri(m.mat, ipiv, work, -1)
if int(work[0]) > len(work) { if int(work[0]) > len(work) {
l := int(work[0]) l := int(work[0])
putFloats(work) putFloat64s(work)
work = getFloats(l, false) work = getFloat64s(l, false)
} }
defer putFloats(work) defer putFloat64s(work)
ok = lapack64.Getri(m.mat, ipiv, work, len(work)) ok = lapack64.Getri(m.mat, ipiv, work, len(work))
if !ok || rcond == 0 { if !ok || rcond == 0 {
// A is exactly singular. // A is exactly singular.
@@ -321,10 +321,10 @@ func (m *Dense) Mul(a, b Matrix) {
case *SymDense: case *SymDense:
if aTrans { if aTrans {
c := getWorkspace(ac, ar, false) c := getDenseWorkspace(ac, ar, false)
blas64.Symm(blas.Left, 1, bU.mat, aU.mat, 0, c.mat) blas64.Symm(blas.Left, 1, bU.mat, aU.mat, 0, c.mat)
strictCopy(m, c.T()) strictCopy(m, c.T())
putWorkspace(c) putDenseWorkspace(c)
return return
} }
blas64.Symm(blas.Right, 1, bU.mat, aU.mat, 0, m.mat) blas64.Symm(blas.Right, 1, bU.mat, aU.mat, 0, m.mat)
@@ -333,7 +333,7 @@ func (m *Dense) Mul(a, b Matrix) {
case *TriDense: case *TriDense:
// Trmm updates in place, so copy aU first. // Trmm updates in place, so copy aU first.
if aTrans { if aTrans {
c := getWorkspace(ac, ar, false) c := getDenseWorkspace(ac, ar, false)
var tmp Dense var tmp Dense
tmp.SetRawMatrix(aU.mat) tmp.SetRawMatrix(aU.mat)
c.Copy(&tmp) c.Copy(&tmp)
@@ -343,7 +343,7 @@ func (m *Dense) Mul(a, b Matrix) {
} }
blas64.Trmm(blas.Left, bT, 1, bU.mat, c.mat) blas64.Trmm(blas.Left, bT, 1, bU.mat, c.mat)
strictCopy(m, c.T()) strictCopy(m, c.T())
putWorkspace(c) putDenseWorkspace(c)
return return
} }
m.Copy(a) m.Copy(a)
@@ -380,10 +380,10 @@ func (m *Dense) Mul(a, b Matrix) {
switch aU := aU.(type) { switch aU := aU.(type) {
case *SymDense: case *SymDense:
if bTrans { if bTrans {
c := getWorkspace(bc, br, false) c := getDenseWorkspace(bc, br, false)
blas64.Symm(blas.Right, 1, aU.mat, bU.mat, 0, c.mat) blas64.Symm(blas.Right, 1, aU.mat, bU.mat, 0, c.mat)
strictCopy(m, c.T()) strictCopy(m, c.T())
putWorkspace(c) putDenseWorkspace(c)
return return
} }
blas64.Symm(blas.Left, 1, aU.mat, bU.mat, 0, m.mat) blas64.Symm(blas.Left, 1, aU.mat, bU.mat, 0, m.mat)
@@ -392,7 +392,7 @@ func (m *Dense) Mul(a, b Matrix) {
case *TriDense: case *TriDense:
// Trmm updates in place, so copy bU first. // Trmm updates in place, so copy bU first.
if bTrans { if bTrans {
c := getWorkspace(bc, br, false) c := getDenseWorkspace(bc, br, false)
var tmp Dense var tmp Dense
tmp.SetRawMatrix(bU.mat) tmp.SetRawMatrix(bU.mat)
c.Copy(&tmp) c.Copy(&tmp)
@@ -402,7 +402,7 @@ func (m *Dense) Mul(a, b Matrix) {
} }
blas64.Trmm(blas.Right, aT, 1, aU.mat, c.mat) blas64.Trmm(blas.Right, aT, 1, aU.mat, c.mat)
strictCopy(m, c.T()) strictCopy(m, c.T())
putWorkspace(c) putDenseWorkspace(c)
return return
} }
m.Copy(b) m.Copy(b)
@@ -441,8 +441,8 @@ func (m *Dense) Mul(a, b Matrix) {
m.checkOverlapMatrix(aU) m.checkOverlapMatrix(aU)
m.checkOverlapMatrix(bU) m.checkOverlapMatrix(bU)
row := getFloats(ac, false) row := getFloat64s(ac, false)
defer putFloats(row) defer putFloat64s(row)
for r := 0; r < ar; r++ { for r := 0; r < ar; r++ {
for i := range row { for i := range row {
row[i] = a.At(r, i) row[i] = a.At(r, i)
@@ -504,19 +504,19 @@ func (m *Dense) Exp(a Matrix) {
a1 := m a1 := m
a1.Copy(a) a1.Copy(a)
v := getWorkspace(r, r, true) v := getDenseWorkspace(r, r, true)
vraw := v.RawMatrix() vraw := v.RawMatrix()
n := r * r n := r * r
vvec := blas64.Vector{N: n, Inc: 1, Data: vraw.Data} vvec := blas64.Vector{N: n, Inc: 1, Data: vraw.Data}
defer putWorkspace(v) defer putDenseWorkspace(v)
u := getWorkspace(r, r, true) u := getDenseWorkspace(r, r, true)
uraw := u.RawMatrix() uraw := u.RawMatrix()
uvec := blas64.Vector{N: n, Inc: 1, Data: uraw.Data} uvec := blas64.Vector{N: n, Inc: 1, Data: uraw.Data}
defer putWorkspace(u) defer putDenseWorkspace(u)
a2 := getWorkspace(r, r, false) a2 := getDenseWorkspace(r, r, false)
defer putWorkspace(a2) defer putDenseWorkspace(a2)
n1 := Norm(a, 1) n1 := Norm(a, 1)
for i, t := range pade { for i, t := range pade {
@@ -526,10 +526,10 @@ func (m *Dense) Exp(a Matrix) {
// This loop only executes once, so // This loop only executes once, so
// this is not as horrible as it looks. // this is not as horrible as it looks.
p := getWorkspace(r, r, true) p := getDenseWorkspace(r, r, true)
praw := p.RawMatrix() praw := p.RawMatrix()
pvec := blas64.Vector{N: n, Inc: 1, Data: praw.Data} pvec := blas64.Vector{N: n, Inc: 1, Data: praw.Data}
defer putWorkspace(p) defer putDenseWorkspace(p)
for k := 0; k < r; k++ { for k := 0; k < r; k++ {
p.set(k, k, 1) p.set(k, k, 1)
@@ -571,27 +571,27 @@ func (m *Dense) Exp(a Matrix) {
} }
a2.Mul(a1, a1) a2.Mul(a1, a1)
i := getWorkspace(r, r, true) i := getDenseWorkspace(r, r, true)
for j := 0; j < r; j++ { for j := 0; j < r; j++ {
i.set(j, j, 1) i.set(j, j, 1)
} }
iraw := i.RawMatrix() iraw := i.RawMatrix()
ivec := blas64.Vector{N: n, Inc: 1, Data: iraw.Data} ivec := blas64.Vector{N: n, Inc: 1, Data: iraw.Data}
defer putWorkspace(i) defer putDenseWorkspace(i)
a2raw := a2.RawMatrix() a2raw := a2.RawMatrix()
a2vec := blas64.Vector{N: n, Inc: 1, Data: a2raw.Data} a2vec := blas64.Vector{N: n, Inc: 1, Data: a2raw.Data}
a4 := getWorkspace(r, r, false) a4 := getDenseWorkspace(r, r, false)
a4raw := a4.RawMatrix() a4raw := a4.RawMatrix()
a4vec := blas64.Vector{N: n, Inc: 1, Data: a4raw.Data} a4vec := blas64.Vector{N: n, Inc: 1, Data: a4raw.Data}
defer putWorkspace(a4) defer putDenseWorkspace(a4)
a4.Mul(a2, a2) a4.Mul(a2, a2)
a6 := getWorkspace(r, r, false) a6 := getDenseWorkspace(r, r, false)
a6raw := a6.RawMatrix() a6raw := a6.RawMatrix()
a6vec := blas64.Vector{N: n, Inc: 1, Data: a6raw.Data} a6vec := blas64.Vector{N: n, Inc: 1, Data: a6raw.Data}
defer putWorkspace(a6) defer putDenseWorkspace(a6)
a6.Mul(a2, a4) a6.Mul(a2, a4)
// V = A_6(b_12*A_6 + b_10*A_4 + b_8*A_2) + b_6*A_6 + b_4*A_4 + b_2*A_2 +b_0*I // V = A_6(b_12*A_6 + b_10*A_4 + b_8*A_2) + b_6*A_6 + b_4*A_4 + b_2*A_2 +b_0*I
@@ -659,11 +659,11 @@ func (m *Dense) Pow(a Matrix, n int) {
} }
// Perform iterative exponentiation by squaring in work space. // Perform iterative exponentiation by squaring in work space.
w := getWorkspace(r, r, false) w := getDenseWorkspace(r, r, false)
w.Copy(a) w.Copy(a)
s := getWorkspace(r, r, false) s := getDenseWorkspace(r, r, false)
s.Copy(a) s.Copy(a)
x := getWorkspace(r, r, false) x := getDenseWorkspace(r, r, false)
for n--; n > 0; n >>= 1 { for n--; n > 0; n >>= 1 {
if n&1 != 0 { if n&1 != 0 {
x.Mul(w, s) x.Mul(w, s)
@@ -675,9 +675,9 @@ func (m *Dense) Pow(a Matrix, n int) {
} }
} }
m.Copy(w) m.Copy(w)
putWorkspace(w) putDenseWorkspace(w)
putWorkspace(s) putDenseWorkspace(s)
putWorkspace(x) putDenseWorkspace(x)
} }
// Kronecker calculates the Kronecker product of a and b, placing the result in // Kronecker calculates the Kronecker product of a and b, placing the result in

View File

@@ -50,9 +50,9 @@ func (e *EigenSym) Factorize(a Symmetric, vectors bool) (ok bool) {
work := []float64{0} work := []float64{0}
lapack64.Syev(jobz, sd.mat, w, work, -1) lapack64.Syev(jobz, sd.mat, w, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
ok = lapack64.Syev(jobz, sd.mat, w, work, len(work)) ok = lapack64.Syev(jobz, sd.mat, w, work, len(work))
putFloats(work) putFloat64s(work)
if !ok { if !ok {
e.vectorsComputed = false e.vectorsComputed = false
e.values = nil e.values = nil
@@ -195,16 +195,16 @@ func (e *Eigen) Factorize(a Matrix, kind EigenKind) (ok bool) {
jobvr = lapack.RightEVCompute jobvr = lapack.RightEVCompute
} }
wr := getFloats(c, false) wr := getFloat64s(c, false)
defer putFloats(wr) defer putFloat64s(wr)
wi := getFloats(c, false) wi := getFloat64s(c, false)
defer putFloats(wi) defer putFloat64s(wi)
work := []float64{0} work := []float64{0}
lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, -1) lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
first := lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, len(work)) first := lapack64.Geev(jobvl, jobvr, sd.mat, wr, wi, vl.mat, vr.mat, work, len(work))
putFloats(work) putFloat64s(work)
if first != 0 { if first != 0 {
e.values = nil e.values = nil

View File

@@ -79,10 +79,10 @@ func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) {
} }
} }
s := getWorkspace(c, c, true) s := getDenseWorkspace(c, c, true)
defer putWorkspace(s) defer putDenseWorkspace(s)
sij := getWorkspace(c, c, false) sij := getDenseWorkspace(c, c, false)
defer putWorkspace(sij) defer putDenseWorkspace(sij)
for i, ai := range a { for i, ai := range a {
for _, aj := range a[i+1:] { for _, aj := range a[i+1:] {
gsvd.err = ai.SolveCholTo(sij, &aj) gsvd.err = ai.SolveCholTo(sij, &aj)
@@ -126,8 +126,8 @@ func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) {
} }
b := make([]Dense, len(m)) b := make([]Dense, len(m))
biT := getWorkspace(c, r, false) biT := getDenseWorkspace(c, r, false)
defer putWorkspace(biT) defer putDenseWorkspace(biT)
for i, d := range m { for i, d := range m {
// All calls to reset will leave an emptied // All calls to reset will leave an emptied
// matrix with capacity to store the result // matrix with capacity to store the result

View File

@@ -31,12 +31,12 @@ func (lq *LQ) updateCond(norm lapack.MatrixNorm) {
// is not the case for CondNorm. Hopefully the error is negligible: κ // is not the case for CondNorm. Hopefully the error is negligible: κ
// is only a qualitative measure anyway. // is only a qualitative measure anyway.
m := lq.lq.mat.Rows m := lq.lq.mat.Rows
work := getFloats(3*m, false) work := getFloat64s(3*m, false)
iwork := getInts(m, false) iwork := getInts(m, false)
l := lq.lq.asTriDense(m, blas.NonUnit, blas.Lower) l := lq.lq.asTriDense(m, blas.NonUnit, blas.Lower)
v := lapack64.Trcon(norm, l.mat, work, iwork) v := lapack64.Trcon(norm, l.mat, work, iwork)
lq.cond = 1 / v lq.cond = 1 / v
putFloats(work) putFloat64s(work)
putInts(iwork) putInts(iwork)
} }
@@ -63,9 +63,9 @@ func (lq *LQ) factorize(a Matrix, norm lapack.MatrixNorm) {
work := []float64{0} work := []float64{0}
lq.tau = make([]float64, k) lq.tau = make([]float64, k)
lapack64.Gelqf(lq.lq.mat, lq.tau, work, -1) lapack64.Gelqf(lq.lq.mat, lq.tau, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Gelqf(lq.lq.mat, lq.tau, work, len(work)) lapack64.Gelqf(lq.lq.mat, lq.tau, work, len(work))
putFloats(work) putFloat64s(work)
lq.updateCond(norm) lq.updateCond(norm)
} }
@@ -159,9 +159,9 @@ func (lq *LQ) QTo(dst *Dense) {
// Construct Q from the elementary reflectors. // Construct Q from the elementary reflectors.
work := []float64{0} work := []float64{0}
lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, -1) lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, len(work)) lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, len(work))
putFloats(work) putFloat64s(work)
} }
// SolveTo finds a minimum-norm solution to a system of linear equations defined // SolveTo finds a minimum-norm solution to a system of linear equations defined
@@ -199,15 +199,15 @@ func (lq *LQ) SolveTo(dst *Dense, trans bool, b Matrix) error {
} }
// Do not need to worry about overlap between x and b because w has its own // Do not need to worry about overlap between x and b because w has its own
// independent storage. // independent storage.
w := getWorkspace(max(r, c), bc, false) w := getDenseWorkspace(max(r, c), bc, false)
w.Copy(b) w.Copy(b)
t := lq.lq.asTriDense(lq.lq.mat.Rows, blas.NonUnit, blas.Lower).mat t := lq.lq.asTriDense(lq.lq.mat.Rows, blas.NonUnit, blas.Lower).mat
if trans { if trans {
work := []float64{0} work := []float64{0}
lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, w.mat, work, -1) lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, w.mat, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, w.mat, work, len(work)) lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, w.mat, work, len(work))
putFloats(work) putFloat64s(work)
ok := lapack64.Trtrs(blas.Trans, t, w.mat) ok := lapack64.Trtrs(blas.Trans, t, w.mat)
if !ok { if !ok {
@@ -223,13 +223,13 @@ func (lq *LQ) SolveTo(dst *Dense, trans bool, b Matrix) error {
} }
work := []float64{0} work := []float64{0}
lapack64.Ormlq(blas.Left, blas.Trans, lq.lq.mat, lq.tau, w.mat, work, -1) lapack64.Ormlq(blas.Left, blas.Trans, lq.lq.mat, lq.tau, w.mat, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Ormlq(blas.Left, blas.Trans, lq.lq.mat, lq.tau, w.mat, work, len(work)) lapack64.Ormlq(blas.Left, blas.Trans, lq.lq.mat, lq.tau, w.mat, work, len(work))
putFloats(work) putFloat64s(work)
} }
// x was set above to be the correct size for the result. // x was set above to be the correct size for the result.
dst.Copy(w) dst.Copy(w)
putWorkspace(w) putDenseWorkspace(w)
if lq.cond > ConditionTolerance { if lq.cond > ConditionTolerance {
return Condition(lq.cond) return Condition(lq.cond)
} }

View File

@@ -30,8 +30,8 @@ type LU struct {
// norm of the original matrix. If anorm is negative it will be estimated. // norm of the original matrix. If anorm is negative it will be estimated.
func (lu *LU) updateCond(anorm float64, norm lapack.MatrixNorm) { func (lu *LU) updateCond(anorm float64, norm lapack.MatrixNorm) {
n := lu.lu.mat.Cols n := lu.lu.mat.Cols
work := getFloats(4*n, false) work := getFloat64s(4*n, false)
defer putFloats(work) defer putFloat64s(work)
iwork := getInts(n, false) iwork := getInts(n, false)
defer putInts(iwork) defer putInts(iwork)
if anorm < 0 { if anorm < 0 {
@@ -79,9 +79,9 @@ func (lu *LU) factorize(a Matrix, norm lapack.MatrixNorm) {
lu.pivot = make([]int, r) lu.pivot = make([]int, r)
} }
lu.pivot = lu.pivot[:r] lu.pivot = lu.pivot[:r]
work := getFloats(r, false) work := getFloat64s(r, false)
anorm := lapack64.Lange(norm, lu.lu.mat, work) anorm := lapack64.Lange(norm, lu.lu.mat, work)
putFloats(work) putFloat64s(work)
lapack64.Getrf(lu.lu.mat, lu.pivot) lapack64.Getrf(lu.lu.mat, lu.pivot)
lu.updateCond(anorm, norm) lu.updateCond(anorm, norm)
} }
@@ -131,8 +131,8 @@ func (lu *LU) LogDet() (det float64, sign float64) {
} }
_, n := lu.lu.Dims() _, n := lu.lu.Dims()
logDiag := getFloats(n, false) logDiag := getFloat64s(n, false)
defer putFloats(logDiag) defer putFloat64s(logDiag)
sign = 1.0 sign = 1.0
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
v := lu.lu.at(i, i) v := lu.lu.at(i, i)
@@ -214,10 +214,10 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y Vector) {
lu.lu.Copy(orig.lu) lu.lu.Copy(orig.lu)
} }
xs := getFloats(n, false) xs := getFloat64s(n, false)
defer putFloats(xs) defer putFloat64s(xs)
ys := getFloats(n, false) ys := getFloat64s(n, false)
defer putFloats(ys) defer putFloat64s(ys)
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
xs[i] = x.AtVec(i) xs[i] = x.AtVec(i)
ys[i] = y.AtVec(i) ys[i] = y.AtVec(i)

View File

@@ -94,10 +94,10 @@ func init() {
} }
} }
// getWorkspace returns a *Dense of size r×c and a data slice // getDenseWorkspace returns a *Dense of size r×c and a data slice
// with a cap that is less than 2*r*c. If clear is true, the // with a cap that is less than 2*r*c. If clear is true, the
// data slice visible through the Matrix interface is zeroed. // data slice visible through the Matrix interface is zeroed.
func getWorkspace(r, c int, clear bool) *Dense { func getDenseWorkspace(r, c int, clear bool) *Dense {
l := uint(r * c) l := uint(r * c)
w := poolDense[poolFor(l)].Get().(*Dense) w := poolDense[poolFor(l)].Get().(*Dense)
w.mat.Data = w.mat.Data[:l] w.mat.Data = w.mat.Data[:l]
@@ -112,17 +112,17 @@ func getWorkspace(r, c int, clear bool) *Dense {
return w return w
} }
// putWorkspace replaces a used *Dense into the appropriate size // putDenseWorkspace replaces a used *Dense into the appropriate size
// workspace pool. putWorkspace must not be called with a matrix // workspace pool. putDenseWorkspace must not be called with a matrix
// where references to the underlying data slice have been kept. // where references to the underlying data slice have been kept.
func putWorkspace(w *Dense) { func putDenseWorkspace(w *Dense) {
poolDense[poolFor(uint(cap(w.mat.Data)))].Put(w) poolDense[poolFor(uint(cap(w.mat.Data)))].Put(w)
} }
// getWorkspaceSym returns a *SymDense of size n and a cap that // getSymDenseWorkspace returns a *SymDense of size n and a cap that
// is less than 2*n. If clear is true, the data slice visible // is less than 2*n. If clear is true, the data slice visible
// through the Matrix interface is zeroed. // through the Matrix interface is zeroed.
func getWorkspaceSym(n int, clear bool) *SymDense { func getSymDenseWorkspace(n int, clear bool) *SymDense {
l := uint(n) l := uint(n)
l *= l l *= l
s := poolSymDense[poolFor(l)].Get().(*SymDense) s := poolSymDense[poolFor(l)].Get().(*SymDense)
@@ -136,17 +136,17 @@ func getWorkspaceSym(n int, clear bool) *SymDense {
return s return s
} }
// putWorkspaceSym replaces a used *SymDense into the appropriate size // putSymDenseWorkspace replaces a used *SymDense into the appropriate size
// workspace pool. putWorkspaceSym must not be called with a matrix // workspace pool. putSymDenseWorkspace must not be called with a matrix
// where references to the underlying data slice have been kept. // where references to the underlying data slice have been kept.
func putWorkspaceSym(s *SymDense) { func putSymDenseWorkspace(s *SymDense) {
poolSymDense[poolFor(uint(cap(s.mat.Data)))].Put(s) poolSymDense[poolFor(uint(cap(s.mat.Data)))].Put(s)
} }
// getWorkspaceTri returns a *TriDense of size n and a cap that // getTriDenseWorkspace returns a *TriDense of size n and a cap that
// is less than 2*n. If clear is true, the data slice visible // is less than 2*n. If clear is true, the data slice visible
// through the Matrix interface is zeroed. // through the Matrix interface is zeroed.
func getWorkspaceTri(n int, kind TriKind, clear bool) *TriDense { func getTriDenseWorkspace(n int, kind TriKind, clear bool) *TriDense {
l := uint(n) l := uint(n)
l *= l l *= l
t := poolTriDense[poolFor(l)].Get().(*TriDense) t := poolTriDense[poolFor(l)].Get().(*TriDense)
@@ -168,17 +168,17 @@ func getWorkspaceTri(n int, kind TriKind, clear bool) *TriDense {
return t return t
} }
// putWorkspaceTri replaces a used *TriDense into the appropriate size // putTriWorkspace replaces a used *TriDense into the appropriate size
// workspace pool. putWorkspaceTri must not be called with a matrix // workspace pool. putTriWorkspace must not be called with a matrix
// where references to the underlying data slice have been kept. // where references to the underlying data slice have been kept.
func putWorkspaceTri(t *TriDense) { func putTriWorkspace(t *TriDense) {
poolTriDense[poolFor(uint(cap(t.mat.Data)))].Put(t) poolTriDense[poolFor(uint(cap(t.mat.Data)))].Put(t)
} }
// getWorkspaceVec returns a *VecDense of length n and a cap that // getVecDenseWorkspace returns a *VecDense of length n and a cap that
// is less than 2*n. If clear is true, the data slice visible // is less than 2*n. If clear is true, the data slice visible
// through the Matrix interface is zeroed. // through the Matrix interface is zeroed.
func getWorkspaceVec(n int, clear bool) *VecDense { func getVecDenseWorkspace(n int, clear bool) *VecDense {
l := uint(n) l := uint(n)
v := poolVecDense[poolFor(l)].Get().(*VecDense) v := poolVecDense[poolFor(l)].Get().(*VecDense)
v.mat.Data = v.mat.Data[:l] v.mat.Data = v.mat.Data[:l]
@@ -189,16 +189,41 @@ func getWorkspaceVec(n int, clear bool) *VecDense {
return v return v
} }
// putWorkspaceVec replaces a used *VecDense into the appropriate size // putVecDenseWorkspace replaces a used *VecDense into the appropriate size
// workspace pool. putWorkspaceVec must not be called with a matrix // workspace pool. putVecDenseWorkspace must not be called with a matrix
// where references to the underlying data slice have been kept. // where references to the underlying data slice have been kept.
func putWorkspaceVec(v *VecDense) { func putVecDenseWorkspace(v *VecDense) {
poolVecDense[poolFor(uint(cap(v.mat.Data)))].Put(v) poolVecDense[poolFor(uint(cap(v.mat.Data)))].Put(v)
} }
// getFloats returns a []float64 of length l and a cap that is // getCDenseWorkspace returns a *CDense of size r×c and a data slice
// with a cap that is less than 2*r*c. If clear is true, the
// data slice visible through the CMatrix interface is zeroed.
func getCDenseWorkspace(r, c int, clear bool) *CDense {
l := uint(r * c)
w := poolCDense[poolFor(l)].Get().(*CDense)
w.mat.Data = w.mat.Data[:l]
if clear {
zeroC(w.mat.Data)
}
w.mat.Rows = r
w.mat.Cols = c
w.mat.Stride = c
w.capRows = r
w.capCols = c
return w
}
// putCDenseWorkspace replaces a used *CDense into the appropriate size
// workspace pool. putWorkspace must not be called with a matrix
// where references to the underlying data slice have been kept.
func putCDenseWorkspace(w *CDense) {
poolCDense[poolFor(uint(cap(w.mat.Data)))].Put(w)
}
// getFloat64s returns a []float64 of length l and a cap that is
// less than 2*l. If clear is true, the slice visible is zeroed. // less than 2*l. If clear is true, the slice visible is zeroed.
func getFloats(l int, clear bool) []float64 { func getFloat64s(l int, clear bool) []float64 {
w := *poolFloat64s[poolFor(uint(l))].Get().(*[]float64) w := *poolFloat64s[poolFor(uint(l))].Get().(*[]float64)
w = w[:l] w = w[:l]
if clear { if clear {
@@ -207,14 +232,14 @@ func getFloats(l int, clear bool) []float64 {
return w return w
} }
// putFloats replaces a used []float64 into the appropriate size // putFloat64s replaces a used []float64 into the appropriate size
// workspace pool. putFloats must not be called with a slice // workspace pool. putFloat64s must not be called with a slice
// where references to the underlying data have been kept. // where references to the underlying data have been kept.
func putFloats(w []float64) { func putFloat64s(w []float64) {
poolFloat64s[poolFor(uint(cap(w)))].Put(&w) poolFloat64s[poolFor(uint(cap(w)))].Put(&w)
} }
// getInts returns a []ints of length l and a cap that is // getInts returns a []int of length l and a cap that is
// less than 2*l. If clear is true, the slice visible is zeroed. // less than 2*l. If clear is true, the slice visible is zeroed.
func getInts(l int, clear bool) []int { func getInts(l int, clear bool) []int {
w := *poolInts[poolFor(uint(l))].Get().(*[]int) w := *poolInts[poolFor(uint(l))].Get().(*[]int)
@@ -233,28 +258,3 @@ func getInts(l int, clear bool) []int {
func putInts(w []int) { func putInts(w []int) {
poolInts[poolFor(uint(cap(w)))].Put(&w) poolInts[poolFor(uint(cap(w)))].Put(&w)
} }
// getWorkspaceCmplx returns a *CDense of size r×c and a data slice
// with a cap that is less than 2*r*c. If clear is true, the
// data slice visible through the CMatrix interface is zeroed.
func getWorkspaceCmplx(r, c int, clear bool) *CDense {
l := uint(r * c)
w := poolCDense[poolFor(l)].Get().(*CDense)
w.mat.Data = w.mat.Data[:l]
if clear {
zeroC(w.mat.Data)
}
w.mat.Rows = r
w.mat.Cols = c
w.mat.Stride = c
w.capRows = r
w.capCols = c
return w
}
// putWorkspaceCmplx replaces a used *CDense into the appropriate size
// workspace pool. putWorkspace must not be called with a matrix
// where references to the underlying data slice have been kept.
func putWorkspaceCmplx(w *CDense) {
poolCDense[poolFor(uint(cap(w.mat.Data)))].Put(w)
}

View File

@@ -20,7 +20,7 @@ func TestPool(t *testing.T) {
for k := 0; k < 5; k++ { for k := 0; k < 5; k++ {
work := make([]*Dense, rand.Intn(10)+1) work := make([]*Dense, rand.Intn(10)+1)
for l := range work { for l := range work {
w := getWorkspace(i, j, true) w := getDenseWorkspace(i, j, true)
if !reflect.DeepEqual(w.mat, m.mat) { if !reflect.DeepEqual(w.mat, m.mat) {
t.Error("unexpected non-zeroed matrix returned by getWorkspace") t.Error("unexpected non-zeroed matrix returned by getWorkspace")
} }
@@ -37,7 +37,7 @@ func TestPool(t *testing.T) {
work[l] = w work[l] = w
} }
for _, w := range work { for _, w := range work {
putWorkspace(w) putDenseWorkspace(w)
} }
} }
} }
@@ -48,8 +48,8 @@ var benchmat *Dense
func poolBenchmark(n, r, c int, clear bool) { func poolBenchmark(n, r, c int, clear bool) {
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
benchmat = getWorkspace(r, c, clear) benchmat = getDenseWorkspace(r, c, clear)
putWorkspace(benchmat) putDenseWorkspace(benchmat)
} }
} }

View File

@@ -45,7 +45,7 @@ func (m *Dense) Product(factors ...Matrix) {
result := p.multiply() result := p.multiply()
m.reuseAsNonZeroed(result.Dims()) m.reuseAsNonZeroed(result.Dims())
m.Copy(result) m.Copy(result)
putWorkspace(result) putDenseWorkspace(result)
} }
// debugProductWalk enables debugging output for Product. // debugProductWalk enables debugging output for Product.
@@ -154,13 +154,13 @@ func (p *multiplier) multiplySubchain(i, j int) (m Matrix, intermediate bool) {
i, ar, ac, result(aTmp), j, br, bc, result(bTmp)) i, ar, ac, result(aTmp), j, br, bc, result(bTmp))
} }
r := getWorkspace(ar, bc, false) r := getDenseWorkspace(ar, bc, false)
r.Mul(a, b) r.Mul(a, b)
if aTmp { if aTmp {
putWorkspace(a.(*Dense)) putDenseWorkspace(a.(*Dense))
} }
if bTmp { if bTmp {
putWorkspace(b.(*Dense)) putDenseWorkspace(b.(*Dense))
} }
return r, true return r, true
} }

View File

@@ -31,11 +31,11 @@ func (qr *QR) updateCond(norm lapack.MatrixNorm) {
// is not the case for CondNorm. Hopefully the error is negligible: κ // is not the case for CondNorm. Hopefully the error is negligible: κ
// is only a qualitative measure anyway. // is only a qualitative measure anyway.
n := qr.qr.mat.Cols n := qr.qr.mat.Cols
work := getFloats(3*n, false) work := getFloat64s(3*n, false)
iwork := getInts(n, false) iwork := getInts(n, false)
r := qr.qr.asTriDense(n, blas.NonUnit, blas.Upper) r := qr.qr.asTriDense(n, blas.NonUnit, blas.Upper)
v := lapack64.Trcon(norm, r.mat, work, iwork) v := lapack64.Trcon(norm, r.mat, work, iwork)
putFloats(work) putFloat64s(work)
putInts(iwork) putInts(iwork)
qr.cond = 1 / v qr.cond = 1 / v
} }
@@ -63,9 +63,9 @@ func (qr *QR) factorize(a Matrix, norm lapack.MatrixNorm) {
work := []float64{0} work := []float64{0}
qr.tau = make([]float64, k) qr.tau = make([]float64, k)
lapack64.Geqrf(qr.qr.mat, qr.tau, work, -1) lapack64.Geqrf(qr.qr.mat, qr.tau, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Geqrf(qr.qr.mat, qr.tau, work, len(work)) lapack64.Geqrf(qr.qr.mat, qr.tau, work, len(work))
putFloats(work) putFloat64s(work)
qr.updateCond(norm) qr.updateCond(norm)
} }
@@ -154,9 +154,9 @@ func (qr *QR) QTo(dst *Dense) {
// Construct Q from the elementary reflectors. // Construct Q from the elementary reflectors.
work := []float64{0} work := []float64{0}
lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, -1) lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, len(work)) lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, len(work))
putFloats(work) putFloat64s(work)
} }
// SolveTo finds a minimum-norm solution to a system of linear equations defined // SolveTo finds a minimum-norm solution to a system of linear equations defined
@@ -194,7 +194,7 @@ func (qr *QR) SolveTo(dst *Dense, trans bool, b Matrix) error {
} }
// Do not need to worry about overlap between m and b because x has its own // Do not need to worry about overlap between m and b because x has its own
// independent storage. // independent storage.
w := getWorkspace(max(r, c), bc, false) w := getDenseWorkspace(max(r, c), bc, false)
w.Copy(b) w.Copy(b)
t := qr.qr.asTriDense(qr.qr.mat.Cols, blas.NonUnit, blas.Upper).mat t := qr.qr.asTriDense(qr.qr.mat.Cols, blas.NonUnit, blas.Upper).mat
if trans { if trans {
@@ -207,15 +207,15 @@ func (qr *QR) SolveTo(dst *Dense, trans bool, b Matrix) error {
} }
work := []float64{0} work := []float64{0}
lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, w.mat, work, -1) lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, w.mat, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, w.mat, work, len(work)) lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, w.mat, work, len(work))
putFloats(work) putFloat64s(work)
} else { } else {
work := []float64{0} work := []float64{0}
lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, w.mat, work, -1) lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, w.mat, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, w.mat, work, len(work)) lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, w.mat, work, len(work))
putFloats(work) putFloat64s(work)
ok := lapack64.Trtrs(blas.NoTrans, t, w.mat) ok := lapack64.Trtrs(blas.NoTrans, t, w.mat)
if !ok { if !ok {
@@ -224,7 +224,7 @@ func (qr *QR) SolveTo(dst *Dense, trans bool, b Matrix) error {
} }
// X was set above to be the correct size for the result. // X was set above to be the correct size for the result.
dst.Copy(w) dst.Copy(w)
putWorkspace(w) putDenseWorkspace(w)
if qr.cond > ConditionTolerance { if qr.cond > ConditionTolerance {
return Condition(qr.cond) return Condition(qr.cond)
} }

View File

@@ -52,10 +52,10 @@ func (m *Dense) Solve(a, b Matrix) error {
case RawMatrixer: case RawMatrixer:
if m != bU || bTrans { if m != bU || bTrans {
if m == bU || m.checkOverlap(rm.RawMatrix()) { if m == bU || m.checkOverlap(rm.RawMatrix()) {
tmp := getWorkspace(br, bc, false) tmp := getDenseWorkspace(br, bc, false)
tmp.Copy(b) tmp.Copy(b)
m.Copy(tmp) m.Copy(tmp)
putWorkspace(tmp) putDenseWorkspace(tmp)
break break
} }
m.Copy(b) m.Copy(b)
@@ -65,19 +65,19 @@ func (m *Dense) Solve(a, b Matrix) error {
m.Copy(b) m.Copy(b)
} else if bTrans { } else if bTrans {
// m and b share data so Copy cannot be used directly. // m and b share data so Copy cannot be used directly.
tmp := getWorkspace(br, bc, false) tmp := getDenseWorkspace(br, bc, false)
tmp.Copy(b) tmp.Copy(b)
m.Copy(tmp) m.Copy(tmp)
putWorkspace(tmp) putDenseWorkspace(tmp)
} }
} }
rm := rma.RawTriangular() rm := rma.RawTriangular()
blas64.Trsm(side, tA, 1, rm, m.mat) blas64.Trsm(side, tA, 1, rm, m.mat)
work := getFloats(3*rm.N, false) work := getFloat64s(3*rm.N, false)
iwork := getInts(rm.N, false) iwork := getInts(rm.N, false)
cond := lapack64.Trcon(CondNorm, rm, work, iwork) cond := lapack64.Trcon(CondNorm, rm, work, iwork)
putFloats(work) putFloat64s(work)
putInts(iwork) putInts(iwork)
if cond > ConditionTolerance { if cond > ConditionTolerance {
return Condition(cond) return Condition(cond)

View File

@@ -131,9 +131,9 @@ func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) {
work := []float64{0} work := []float64{0}
lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, -1) lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, -1)
work = getFloats(int(work[0]), false) work = getFloat64s(int(work[0]), false)
ok = lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, len(work)) ok = lapack64.Gesvd(jobU, jobVT, aCopy.mat, svd.u, svd.vt, svd.s, work, len(work))
putFloats(work) putFloat64s(work)
if !ok { if !ok {
svd.kind = 0 svd.kind = 0
} }
@@ -318,12 +318,12 @@ func (svd *SVD) SolveTo(dst *Dense, b Matrix, rank int) []float64 {
s := svd.s[:rank] s := svd.s[:rank]
_, bc := b.Dims() _, bc := b.Dims()
c := getWorkspace(svd.u.Cols, bc, false) c := getDenseWorkspace(svd.u.Cols, bc, false)
defer putWorkspace(c) defer putDenseWorkspace(c)
c.Mul(u.T(), b) c.Mul(u.T(), b)
y := getWorkspace(rank, bc, false) y := getDenseWorkspace(rank, bc, false)
defer putWorkspace(y) defer putDenseWorkspace(y)
y.DivElem(c.slice(0, rank, 0, bc), repVector{vec: s, cols: bc}) y.DivElem(c.slice(0, rank, 0, bc), repVector{vec: s, cols: bc})
dst.Mul(vt.slice(0, rank, 0, svd.vt.Cols).T(), y) dst.Mul(vt.slice(0, rank, 0, svd.vt.Cols).T(), y)
@@ -395,12 +395,12 @@ func (svd *SVD) SolveVecTo(dst *VecDense, b Vector, rank int) float64 {
} }
s := svd.s[:rank] s := svd.s[:rank]
c := getWorkspaceVec(svd.u.Cols, false) c := getVecDenseWorkspace(svd.u.Cols, false)
defer putWorkspaceVec(c) defer putVecDenseWorkspace(c)
c.MulVec(u.T(), b) c.MulVec(u.T(), b)
y := getWorkspaceVec(rank, false) y := getVecDenseWorkspace(rank, false)
defer putWorkspaceVec(y) defer putVecDenseWorkspace(y)
y.DivElemVec(c.sliceVec(0, rank), NewVecDense(rank, s)) y.DivElemVec(c.sliceVec(0, rank), NewVecDense(rank, s))
dst.MulVec(vt.slice(0, rank, 0, svd.vt.Cols).T(), y) dst.MulVec(vt.slice(0, rank, 0, svd.vt.Cols).T(), y)

View File

@@ -257,8 +257,8 @@ func (s *SymBandDense) Norm(norm float64) float64 {
} }
lnorm := normLapack(norm, false) lnorm := normLapack(norm, false)
if lnorm == lapack.MaxColumnSum || lnorm == lapack.MaxRowSum { if lnorm == lapack.MaxColumnSum || lnorm == lapack.MaxRowSum {
work := getFloats(s.mat.N, false) work := getFloat64s(s.mat.N, false)
defer putFloats(work) defer putFloat64s(work)
return lapack64.Lansb(lnorm, s.mat, work) return lapack64.Lansb(lnorm, s.mat, work)
} }
return lapack64.Lansb(lnorm, s.mat, nil) return lapack64.Lansb(lnorm, s.mat, nil)
@@ -293,15 +293,15 @@ func (s *SymBandDense) MulVecTo(dst *VecDense, _ bool, x Vector) {
dst.checkOverlap(xVec.mat) dst.checkOverlap(xVec.mat)
blas64.Sbmv(1, s.mat, xVec.mat, 0, dst.mat) blas64.Sbmv(1, s.mat, xVec.mat, 0, dst.mat)
} else { } else {
xCopy := getWorkspaceVec(n, false) xCopy := getVecDenseWorkspace(n, false)
xCopy.CloneFromVec(xVec) xCopy.CloneFromVec(xVec)
blas64.Sbmv(1, s.mat, xCopy.mat, 0, dst.mat) blas64.Sbmv(1, s.mat, xCopy.mat, 0, dst.mat)
putWorkspaceVec(xCopy) putVecDenseWorkspace(xCopy)
} }
} else { } else {
xCopy := getWorkspaceVec(n, false) xCopy := getVecDenseWorkspace(n, false)
xCopy.CloneFromVec(x) xCopy.CloneFromVec(x)
blas64.Sbmv(1, s.mat, xCopy.mat, 0, dst.mat) blas64.Sbmv(1, s.mat, xCopy.mat, 0, dst.mat)
putWorkspaceVec(xCopy) putVecDenseWorkspace(xCopy)
} }
} }

View File

@@ -238,10 +238,10 @@ func (s *SymDense) isolatedWorkspace(a Symmetric) (w *SymDense, restore func())
if n == 0 { if n == 0 {
panic(ErrZeroLength) panic(ErrZeroLength)
} }
w = getWorkspaceSym(n, false) w = getSymDenseWorkspace(n, false)
return w, func() { return w, func() {
s.CopySym(w) s.CopySym(w)
putWorkspaceSym(w) putSymDenseWorkspace(w)
} }
} }
@@ -405,10 +405,10 @@ func (s *SymDense) SymOuterK(alpha float64, x Matrix) {
panic(badSymTriangle) panic(badSymTriangle)
case s.mat.N == n: case s.mat.N == n:
if s == x { if s == x {
w := getWorkspaceSym(n, true) w := getSymDenseWorkspace(n, true)
w.SymRankK(w, alpha, x) w.SymRankK(w, alpha, x)
s.CopySym(w) s.CopySym(w)
putWorkspaceSym(w) putSymDenseWorkspace(w)
} else { } else {
switch r := x.(type) { switch r := x.(type) {
case RawMatrixer: case RawMatrixer:
@@ -592,8 +592,8 @@ func (s *SymDense) Norm(norm float64) float64 {
} }
lnorm := normLapack(norm, false) lnorm := normLapack(norm, false)
if lnorm == lapack.MaxColumnSum || lnorm == lapack.MaxRowSum { if lnorm == lapack.MaxColumnSum || lnorm == lapack.MaxRowSum {
work := getFloats(s.mat.N, false) work := getFloat64s(s.mat.N, false)
defer putFloats(work) defer putFloat64s(work)
return lapack64.Lansy(lnorm, s.mat, work) return lapack64.Lansy(lnorm, s.mat, work)
} }
return lapack64.Lansy(lnorm, s.mat, nil) return lapack64.Lansy(lnorm, s.mat, nil)

View File

@@ -366,10 +366,10 @@ func (t *TriDense) isolatedWorkspace(a Triangular) (w *TriDense, restore func())
if n == 0 { if n == 0 {
panic(ErrZeroLength) panic(ErrZeroLength)
} }
w = getWorkspaceTri(n, kind, false) w = getTriDenseWorkspace(n, kind, false)
return w, func() { return w, func() {
t.Copy(w) t.Copy(w)
putWorkspaceTri(w) putTriWorkspace(w)
} }
} }
@@ -459,10 +459,10 @@ func (t *TriDense) InverseTri(a Triangular) error {
n, _ := a.Triangle() n, _ := a.Triangle()
t.reuseAsNonZeroed(a.Triangle()) t.reuseAsNonZeroed(a.Triangle())
t.Copy(a) t.Copy(a)
work := getFloats(3*n, false) work := getFloat64s(3*n, false)
iwork := getInts(n, false) iwork := getInts(n, false)
cond := lapack64.Trcon(CondNorm, t.mat, work, iwork) cond := lapack64.Trcon(CondNorm, t.mat, work, iwork)
putFloats(work) putFloat64s(work)
putInts(iwork) putInts(iwork)
if math.IsInf(cond, 1) { if math.IsInf(cond, 1) {
return Condition(cond) return Condition(cond)
@@ -634,8 +634,8 @@ func (t *TriDense) Norm(norm float64) float64 {
} }
lnorm := normLapack(norm, false) lnorm := normLapack(norm, false)
if lnorm == lapack.MaxColumnSum { if lnorm == lapack.MaxColumnSum {
work := getFloats(t.mat.N, false) work := getFloat64s(t.mat.N, false)
defer putFloats(work) defer putFloat64s(work)
return lapack64.Lantr(lnorm, t.mat, work) return lapack64.Lantr(lnorm, t.mat, work)
} }
return lapack64.Lantr(lnorm, t.mat, nil) return lapack64.Lantr(lnorm, t.mat, nil)

View File

@@ -461,8 +461,8 @@ func (t *TriBandDense) Norm(norm float64) float64 {
} }
lnorm := normLapack(norm, false) lnorm := normLapack(norm, false)
if lnorm == lapack.MaxColumnSum { if lnorm == lapack.MaxColumnSum {
work := getFloats(t.mat.N, false) work := getFloat64s(t.mat.N, false)
defer putFloats(work) defer putFloat64s(work)
return lapack64.Lantb(lnorm, t.mat, work) return lapack64.Lantb(lnorm, t.mat, work)
} }
return lapack64.Lantb(lnorm, t.mat, nil) return lapack64.Lantb(lnorm, t.mat, nil)

View File

@@ -216,10 +216,10 @@ func (a *Tridiag) MulVecTo(dst *VecDense, trans bool, x Vector) {
dst.checkOverlap(xVec.mat) dst.checkOverlap(xVec.mat)
lapack64.Lagtm(t, 1, a.mat, xVec.asGeneral(), 0, dst.asGeneral()) lapack64.Lagtm(t, 1, a.mat, xVec.asGeneral(), 0, dst.asGeneral())
} else { } else {
xCopy := getWorkspaceVec(n, false) xCopy := getVecDenseWorkspace(n, false)
xCopy.CloneFromVec(x) xCopy.CloneFromVec(x)
lapack64.Lagtm(t, 1, a.mat, xCopy.asGeneral(), 0, dst.asGeneral()) lapack64.Lagtm(t, 1, a.mat, xCopy.asGeneral(), 0, dst.asGeneral())
putWorkspaceVec(xCopy) putVecDenseWorkspace(xCopy)
} }
} }

View File

@@ -773,10 +773,10 @@ func (v *VecDense) isolatedWorkspace(a Vector) (n *VecDense, restore func()) {
if l == 0 { if l == 0 {
panic(ErrZeroLength) panic(ErrZeroLength)
} }
n = getWorkspaceVec(l, false) n = getVecDenseWorkspace(l, false)
return n, func() { return n, func() {
v.CopyVec(n) v.CopyVec(n)
putWorkspaceVec(n) putVecDenseWorkspace(n)
} }
} }