package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-11-03 03:13:27 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

blas64: add length field N to Vector

Browse Source 
blas64: add length field N to Vector

Alongside, fix the implementation of mat.VecDense and mat.Diagonal, as well as other changes needed to fix this change.

Fixes #736.
This commit is contained in:
Brendan Tracey
2018-12-10 08:36:04 +00:00
committed by GitHub
parent 6b862e40c9
commit 572d9101fe
32 changed files with 228 additions and 223 deletions
Download Patch File Download Diff File Expand all files Collapse all files

55
mat/dense_arithmetic.go
View File

@@ -504,12 +504,13 @@ func (m *Dense) Exp(a Matrix) {
a1.Copy(a)
v := getWorkspace(r, r, true)
vraw := v.RawMatrix()
vvec := blas64.Vector{Inc: 1, Data: vraw.Data}
n := r * r
vvec := blas64.Vector{N: n, Inc: 1, Data: vraw.Data}
defer putWorkspace(v)
u := getWorkspace(r, r, true)
uraw := u.RawMatrix()
uvec := blas64.Vector{Inc: 1, Data: uraw.Data}
uvec := blas64.Vector{N: n, Inc: 1, Data: uraw.Data}
defer putWorkspace(u)
a2 := getWorkspace(r, r, false)
@@ -525,7 +526,7 @@ func (m *Dense) Exp(a Matrix) {
// this is not as horrible as it looks.
p := getWorkspace(r, r, true)
praw := p.RawMatrix()
pvec := blas64.Vector{Inc: 1, Data: praw.Data}
pvec := blas64.Vector{N: n, Inc: 1, Data: praw.Data}
defer putWorkspace(p)
for k := 0; k < r; k++ {
@@ -537,8 +538,8 @@ func (m *Dense) Exp(a Matrix) {
a2.Mul(a1, a1)
for j := 0; j <= i; j++ {
p.Mul(p, a2)
blas64.Axpy(r*r, t.b[2*j+2], pvec, vvec)
blas64.Axpy(r*r, t.b[2*j+3], pvec, uvec)
blas64.Axpy(t.b[2*j+2], pvec, vvec)
blas64.Axpy(t.b[2*j+3], pvec, uvec)
}
u.Mul(a1, u)
@@ -573,43 +574,43 @@ func (m *Dense) Exp(a Matrix) {
i.set(j, j, 1)
}
iraw := i.RawMatrix()
ivec := blas64.Vector{Inc: 1, Data: iraw.Data}
ivec := blas64.Vector{N: n, Inc: 1, Data: iraw.Data}
defer putWorkspace(i)
a2raw := a2.RawMatrix()
a2vec := blas64.Vector{Inc: 1, Data: a2raw.Data}
a2vec := blas64.Vector{N: n, Inc: 1, Data: a2raw.Data}
a4 := getWorkspace(r, r, false)
a4raw := a4.RawMatrix()
a4vec := blas64.Vector{Inc: 1, Data: a4raw.Data}
a4vec := blas64.Vector{N: n, Inc: 1, Data: a4raw.Data}
defer putWorkspace(a4)
a4.Mul(a2, a2)
a6 := getWorkspace(r, r, false)
a6raw := a6.RawMatrix()
a6vec := blas64.Vector{Inc: 1, Data: a6raw.Data}
a6vec := blas64.Vector{N: n, Inc: 1, Data: a6raw.Data}
defer putWorkspace(a6)
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
blas64.Axpy(r*r, b[12], a6vec, vvec)
blas64.Axpy(r*r, b[10], a4vec, vvec)
blas64.Axpy(r*r, b[8], a2vec, vvec)
blas64.Axpy(b[12], a6vec, vvec)
blas64.Axpy(b[10], a4vec, vvec)
blas64.Axpy(b[8], a2vec, vvec)
v.Mul(v, a6)
blas64.Axpy(r*r, b[6], a6vec, vvec)
blas64.Axpy(r*r, b[4], a4vec, vvec)
blas64.Axpy(r*r, b[2], a2vec, vvec)
blas64.Axpy(r*r, b[0], ivec, vvec)
blas64.Axpy(b[6], a6vec, vvec)
blas64.Axpy(b[4], a4vec, vvec)
blas64.Axpy(b[2], a2vec, vvec)
blas64.Axpy(b[0], ivec, vvec)
// U = A(A_6(b_13*A_6 + b_11*A_4 + b_9*A_2) + b_7*A_6 + b_5*A_4 + b_2*A_3 +b_1*I)
blas64.Axpy(r*r, b[13], a6vec, uvec)
blas64.Axpy(r*r, b[11], a4vec, uvec)
blas64.Axpy(r*r, b[9], a2vec, uvec)
blas64.Axpy(b[13], a6vec, uvec)
blas64.Axpy(b[11], a4vec, uvec)
blas64.Axpy(b[9], a2vec, uvec)
u.Mul(u, a6)
blas64.Axpy(r*r, b[7], a6vec, uvec)
blas64.Axpy(r*r, b[5], a4vec, uvec)
blas64.Axpy(r*r, b[3], a2vec, uvec)
blas64.Axpy(r*r, b[1], ivec, uvec)
blas64.Axpy(b[7], a6vec, uvec)
blas64.Axpy(b[5], a4vec, uvec)
blas64.Axpy(b[3], a2vec, uvec)
blas64.Axpy(b[1], ivec, uvec)
u.Mul(u, a1)
// Use i as a workspace here and
@@ -783,14 +784,14 @@ func (m *Dense) RankOne(a Matrix, alpha float64, x, y Vector) {
xU, _ := untranspose(x)
if rv, ok := xU.(RawVectorer); ok {
xmat = rv.RawVector()
m.checkOverlap((&VecDense{mat: xmat, n: x.Len()}).asGeneral())
m.checkOverlap((&VecDense{mat: xmat}).asGeneral())
} else {
fast = false
}
yU, _ := untranspose(y)
if rv, ok := yU.(RawVectorer); ok {
ymat = rv.RawVector()
m.checkOverlap((&VecDense{mat: ymat, n: y.Len()}).asGeneral())
m.checkOverlap((&VecDense{mat: ymat}).asGeneral())
} else {
fast = false
}
@@ -856,7 +857,7 @@ func (m *Dense) Outer(alpha float64, x, y Vector) {
xU, _ := untranspose(x)
if rv, ok := xU.(RawVectorer); ok {
xmat = rv.RawVector()
m.checkOverlap((&VecDense{mat: xmat, n: x.Len()}).asGeneral())
m.checkOverlap((&VecDense{mat: xmat}).asGeneral())
} else {
fast = false
@@ -864,7 +865,7 @@ func (m *Dense) Outer(alpha float64, x, y Vector) {
yU, _ := untranspose(y)
if rv, ok := yU.(RawVectorer); ok {
ymat = rv.RawVector()
m.checkOverlap((&VecDense{mat: ymat, n: y.Len()}).asGeneral())
m.checkOverlap((&VecDense{mat: ymat}).asGeneral())
} else {
fast = false
}