Skip to content

Commit

Permalink
Merge 47a1f0d into 6b862e4
Browse files Browse the repository at this point in the history
  • Loading branch information
@btracey
btracey committed Dec 9, 2018
2 parents 6b862e4 + 47a1f0d commit 30781f0
Show file tree
Hide file tree
Showing 32 changed files with 228 additions and 223 deletions.
67 changes: 45 additions & 22 deletions blas/blas64/blas64.go
View file
Expand Up @@ -28,8 +28,9 @@ func Implementation() blas.Float64 {

// Vector represents a vector with an associated element increment.
type Vector struct {
Data []float64
N int
Inc int
Data []float64
}

// General represents a matrix using the conventional storage scheme.
Expand Down Expand Up @@ -98,64 +99,80 @@ type SymmetricPacked struct {

// Level 1

const negInc = "blas64: negative vector increment"
const (
negInc = "blas64: negative vector increment"
badLength = "blas64: vector length mismatch"
)

// Dot computes the dot product of the two vectors:
// \sum_i x[i]*y[i].
func Dot(n int, x, y Vector) float64 {
return blas64.Ddot(n, x.Data, x.Inc, y.Data, y.Inc)
func Dot(x, y Vector) float64 {
if x.N != y.N {
panic(badLength)
}
return blas64.Ddot(x.N, x.Data, x.Inc, y.Data, y.Inc)
}

// Nrm2 computes the Euclidean norm of the vector x:
// sqrt(\sum_i x[i]*x[i]).
//
// Nrm2 will panic if the vector increment is negative.
func Nrm2(n int, x Vector) float64 {
func Nrm2(x Vector) float64 {
if x.Inc < 0 {
panic(negInc)
}
return blas64.Dnrm2(n, x.Data, x.Inc)
return blas64.Dnrm2(x.N, x.Data, x.Inc)
}

// Asum computes the sum of the absolute values of the elements of x:
// \sum_i |x[i]|.
//
// Asum will panic if the vector increment is negative.
func Asum(n int, x Vector) float64 {
func Asum(x Vector) float64 {
if x.Inc < 0 {
panic(negInc)
}
return blas64.Dasum(n, x.Data, x.Inc)
return blas64.Dasum(x.N, x.Data, x.Inc)
}

// Iamax returns the index of an element of x with the largest absolute value.
// If there are multiple such indices the earliest is returned.
// Iamax returns -1 if n == 0.
//
// Iamax will panic if the vector increment is negative.
func Iamax(n int, x Vector) int {
func Iamax(x Vector) int {
if x.Inc < 0 {
panic(negInc)
}
return blas64.Idamax(n, x.Data, x.Inc)
return blas64.Idamax(x.N, x.Data, x.Inc)
}

// Swap exchanges the elements of the two vectors:
// x[i], y[i] = y[i], x[i] for all i.
func Swap(n int, x, y Vector) {
blas64.Dswap(n, x.Data, x.Inc, y.Data, y.Inc)
func Swap(x, y Vector) {
if x.N != y.N {
panic(badLength)
}
blas64.Dswap(x.N, x.Data, x.Inc, y.Data, y.Inc)
}

// Copy copies the elements of x into the elements of y:
// y[i] = x[i] for all i.
func Copy(n int, x, y Vector) {
blas64.Dcopy(n, x.Data, x.Inc, y.Data, y.Inc)
// Copy requires that the lengths of x and y match and will panic otherwise.
func Copy(x, y Vector) {
if x.N != y.N {
panic(badLength)
}
blas64.Dcopy(x.N, x.Data, x.Inc, y.Data, y.Inc)
}

// Axpy adds x scaled by alpha to y:
// y[i] += alpha*x[i] for all i.
func Axpy(n int, alpha float64, x, y Vector) {
blas64.Daxpy(n, alpha, x.Data, x.Inc, y.Data, y.Inc)
func Axpy(alpha float64, x, y Vector) {
if x.N != y.N {
panic(badLength)
}
blas64.Daxpy(x.N, alpha, x.Data, x.Inc, y.Data, y.Inc)
}

// Rotg computes the parameters of a Givens plane rotation so that
Expand Down Expand Up @@ -185,25 +202,31 @@ func Rotmg(d1, d2, b1, b2 float64) (p blas.DrotmParams, rd1, rd2, rb1 float64) {
// and y:
// x[i] = c*x[i] + s*y[i],
// y[i] = -s*x[i] + c*y[i], for all i.
func Rot(n int, x, y Vector, c, s float64) {
blas64.Drot(n, x.Data, x.Inc, y.Data, y.Inc, c, s)
func Rot(x, y Vector, c, s float64) {
if x.N != y.N {
panic(badLength)
}
blas64.Drot(x.N, x.Data, x.Inc, y.Data, y.Inc, c, s)
}

// Rotm applies the modified Givens rotation to n points represented by the
// vectors x and y.
func Rotm(n int, x, y Vector, p blas.DrotmParams) {
blas64.Drotm(n, x.Data, x.Inc, y.Data, y.Inc, p)
func Rotm(x, y Vector, p blas.DrotmParams) {
if x.N != y.N {
panic(badLength)
}
blas64.Drotm(x.N, x.Data, x.Inc, y.Data, y.Inc, p)
}

// Scal scales the vector x by alpha:
// x[i] *= alpha for all i.
//
// Scal will panic if the vector increment is negative.
func Scal(n int, alpha float64, x Vector) {
func Scal(alpha float64, x Vector) {
if x.Inc < 0 {
panic(negInc)
}
blas64.Dscal(n, alpha, x.Data, x.Inc)
blas64.Dscal(x.N, alpha, x.Data, x.Inc)
}

// Level 2
Expand Down
10 changes: 4 additions & 6 deletions lapack/testlapack/dgebak.go
View file
Expand Up @@ -72,15 +72,13 @@ func testDgebak(t *testing.T, impl Dgebaker, job lapack.BalanceJob, side lapack.
// Make up some random permutations.
for i := n - 1; i > ihi; i-- {
scale[i] = float64(rnd.Intn(i + 1))
blas64.Swap(n,
blas64.Vector{Data: p.Data[i:], Inc: p.Stride},
blas64.Vector{Data: p.Data[int(scale[i]):], Inc: p.Stride})
blas64.Swap(blas64.Vector{N: n, Data: p.Data[i:], Inc: p.Stride},
blas64.Vector{N: n, Data: p.Data[int(scale[i]):], Inc: p.Stride})
}
for i := 0; i < ilo; i++ {
scale[i] = float64(i + rnd.Intn(ihi-i+1))
blas64.Swap(n,
blas64.Vector{Data: p.Data[i:], Inc: p.Stride},
blas64.Vector{Data: p.Data[int(scale[i]):], Inc: p.Stride})
blas64.Swap(blas64.Vector{N: n, Data: p.Data[i:], Inc: p.Stride},
blas64.Vector{N: n, Data: p.Data[int(scale[i]):], Inc: p.Stride})
}
}

Expand Down
10 changes: 4 additions & 6 deletions lapack/testlapack/dgebal.go
View file
Expand Up @@ -144,14 +144,12 @@ func testDgebal(t *testing.T, impl Dgebaler, job lapack.BalanceJob, a blas64.Gen
// Create the permutation matrix P.
p := eye(n, n)
for j := n - 1; j > ihi; j-- {
blas64.Swap(n,
blas64.Vector{Data: p.Data[j:], Inc: p.Stride},
blas64.Vector{Data: p.Data[int(scale[j]):], Inc: p.Stride})
blas64.Swap(blas64.Vector{N: n, Data: p.Data[j:], Inc: p.Stride},
blas64.Vector{N: n, Data: p.Data[int(scale[j]):], Inc: p.Stride})
}
for j := 0; j < ilo; j++ {
blas64.Swap(n,
blas64.Vector{Data: p.Data[j:], Inc: p.Stride},
blas64.Vector{Data: p.Data[int(scale[j]):], Inc: p.Stride})
blas64.Swap(blas64.Vector{N: n, Data: p.Data[j:], Inc: p.Stride},
blas64.Vector{N: n, Data: p.Data[int(scale[j]):], Inc: p.Stride})
}
// Compute P^T*A*P and store into want.
ap := zeros(n, n, n)
Expand Down
3 changes: 2 additions & 1 deletion lapack/testlapack/dgetf2.go
View file
Expand Up @@ -170,7 +170,8 @@ func checkPLU(t *testing.T, ok bool, m, n, lda int, ipiv []int, factorized, orig
}
for i := len(ipiv) - 1; i >= 0; i-- {
v := ipiv[i]
blas64.Swap(m, blas64.Vector{Inc: 1, Data: p[i*ldp:]}, blas64.Vector{Inc: 1, Data: p[v*ldp:]})
blas64.Swap(blas64.Vector{N: m, Inc: 1, Data: p[i*ldp:]},
blas64.Vector{N: m, Inc: 1, Data: p[v*ldp:]})
}
P := blas64.General{
Rows: m,
Expand Down
8 changes: 4 additions & 4 deletions lapack/testlapack/dlange.go
View file
Expand Up @@ -49,7 +49,7 @@ func DlangeTest(t *testing.T, impl Dlanger) {
norm := impl.Dlange(lapack.MaxAbs, m, n, a, lda, work)
var ans float64
for i := 0; i < m; i++ {
idx := blas64.Iamax(n, blas64.Vector{Inc: 1, Data: aCopy[i*lda:]})
idx := blas64.Iamax(blas64.Vector{N: n, Inc: 1, Data: aCopy[i*lda:]})
ans = math.Max(ans, math.Abs(a[i*lda+idx]))
}
// Should be strictly equal because there is no floating point summation error.
Expand All @@ -61,7 +61,7 @@ func DlangeTest(t *testing.T, impl Dlanger) {
norm = impl.Dlange(lapack.MaxColumnSum, m, n, a, lda, work)
ans = 0
for i := 0; i < n; i++ {
sum := blas64.Asum(m, blas64.Vector{Inc: lda, Data: aCopy[i:]})
sum := blas64.Asum(blas64.Vector{N: m, Inc: lda, Data: aCopy[i:]})
ans = math.Max(ans, sum)
}
if math.Abs(norm-ans) > 1e-14 {
Expand All @@ -72,7 +72,7 @@ func DlangeTest(t *testing.T, impl Dlanger) {
norm = impl.Dlange(lapack.MaxRowSum, m, n, a, lda, work)
ans = 0
for i := 0; i < m; i++ {
sum := blas64.Asum(n, blas64.Vector{Inc: 1, Data: aCopy[i*lda:]})
sum := blas64.Asum(blas64.Vector{N: n, Inc: 1, Data: aCopy[i*lda:]})
ans = math.Max(ans, sum)
}
if math.Abs(norm-ans) > 1e-14 {
Expand All @@ -83,7 +83,7 @@ func DlangeTest(t *testing.T, impl Dlanger) {
norm = impl.Dlange(lapack.Frobenius, m, n, a, lda, work)
ans = 0
for i := 0; i < m; i++ {
sum := blas64.Nrm2(n, blas64.Vector{Inc: 1, Data: aCopy[i*lda:]})
sum := blas64.Nrm2(blas64.Vector{N: n, Inc: 1, Data: aCopy[i*lda:]})
ans += sum * sum
}
ans = math.Sqrt(ans)
Expand Down
16 changes: 6 additions & 10 deletions lapack/testlapack/general.go
View file
Expand Up @@ -828,18 +828,16 @@ func isOrthogonal(q blas64.General) bool {
// an orthonormal basis of the Euclidean space R^n.
const tol = 1e-13
for i := 0; i < n; i++ {
nrm := blas64.Nrm2(n, blas64.Vector{Data: q.Data[i*q.Stride:], Inc: 1})
nrm := blas64.Nrm2(blas64.Vector{N: n, Data: q.Data[i*q.Stride:], Inc: 1})
if math.IsNaN(nrm) {
return false
}
if math.Abs(nrm-1) > tol {
return false
}
for j := i + 1; j < n; j++ {
dot := blas64.Dot(n,
blas64.Vector{Data: q.Data[i*q.Stride:], Inc: 1},
blas64.Vector{Data: q.Data[j*q.Stride:], Inc: 1},
)
dot := blas64.Dot(blas64.Vector{N: n, Data: q.Data[i*q.Stride:], Inc: 1},
blas64.Vector{N: n, Data: q.Data[j*q.Stride:], Inc: 1})
if math.IsNaN(dot) {
return false
}
Expand All @@ -860,18 +858,16 @@ func hasOrthonormalColumns(q blas64.General) bool {
ldq := q.Stride
const tol = 1e-13
for i := 0; i < n; i++ {
nrm := blas64.Nrm2(m, blas64.Vector{Data: q.Data[i:], Inc: ldq})
nrm := blas64.Nrm2(blas64.Vector{N: m, Data: q.Data[i:], Inc: ldq})
if math.IsNaN(nrm) {
return false
}
if math.Abs(nrm-1) > tol {
return false
}
for j := i + 1; j < n; j++ {
dot := blas64.Dot(m,
blas64.Vector{Data: q.Data[i:], Inc: ldq},
blas64.Vector{Data: q.Data[j:], Inc: ldq},
)
dot := blas64.Dot(blas64.Vector{N: m, Data: q.Data[i:], Inc: ldq},
blas64.Vector{N: m, Data: q.Data[j:], Inc: ldq})
if math.IsNaN(dot) {
return false
}
Expand Down
2 changes: 1 addition & 1 deletion mat/band.go
View file
Expand Up @@ -193,10 +193,10 @@ func (b *BandDense) DiagView() Diagonal {
n := min(b.mat.Rows, b.mat.Cols)
return &DiagDense{
mat: blas64.Vector{
N: n,
Inc: b.mat.Stride,
Data: b.mat.Data[b.mat.KL : (n-1)*b.mat.Stride+b.mat.KL+1],
},
n: n,
}
}

Expand Down
22 changes: 11 additions & 11 deletions mat/cholesky.go
View file
Expand Up @@ -478,12 +478,12 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
tmp.CopyVec(x)
xmat = tmp.RawVector()
}
blas64.Copy(n, xmat, blas64.Vector{Data: work, Inc: 1})
blas64.Copy(xmat, blas64.Vector{N: n, Data: work, Inc: 1})

if alpha > 0 {
// Compute rank-1 update.
if alpha != 1 {
blas64.Scal(n, math.Sqrt(alpha), blas64.Vector{Data: work, Inc: 1})
blas64.Scal(math.Sqrt(alpha), blas64.Vector{N: n, Data: work, Inc: 1})
}
umat := c.chol.mat
stride := umat.Stride
Expand All @@ -503,9 +503,9 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
// Multiply the extended factorization matrix by
// the Givens matrix from the left. Only
// the i-th row and x are modified.
blas64.Rot(n-i-1,
blas64.Vector{Data: umat.Data[i*stride+i+1 : i*stride+n], Inc: 1},
blas64.Vector{Data: work[i+1 : n], Inc: 1},
blas64.Rot(
blas64.Vector{N: n - i - 1, Data: umat.Data[i*stride+i+1 : i*stride+n], Inc: 1},
blas64.Vector{N: n - i - 1, Data: work[i+1 : n], Inc: 1},
c, s)
}
}
Expand All @@ -516,7 +516,7 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
// Compute rank-1 downdate.
alpha = math.Sqrt(-alpha)
if alpha != 1 {
blas64.Scal(n, alpha, blas64.Vector{Data: work, Inc: 1})
blas64.Scal(alpha, blas64.Vector{N: n, Data: work, Inc: 1})
}
// Solve U^T * p = x storing the result into work.
ok = lapack64.Trtrs(blas.Trans, c.chol.RawTriangular(), blas64.General{
Expand All @@ -530,7 +530,7 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
// the factorization is valid.
panic(badCholesky)
}
norm := blas64.Nrm2(n, blas64.Vector{Data: work, Inc: 1})
norm := blas64.Nrm2(blas64.Vector{N: n, Data: work, Inc: 1})
if norm >= 1 {
// The updated matrix is not positive definite.
return false
Expand All @@ -557,9 +557,9 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
// Apply Givens matrices to U.
// TODO(vladimir-ch): Use workspace to avoid modifying the
// receiver in case an invalid factorization is created.
blas64.Rot(n-i,
blas64.Vector{Data: work[i:n], Inc: 1},
blas64.Vector{Data: umat.Data[i*stride+i : i*stride+n], Inc: 1},
blas64.Rot(
blas64.Vector{N: n - i, Data: work[i:n], Inc: 1},
blas64.Vector{N: n - i, Data: umat.Data[i*stride+i : i*stride+n], Inc: 1},
cos[i], sin[i])
if umat.Data[i*stride+i] == 0 {
// The matrix is singular (may rarely happen due to
Expand All @@ -570,7 +570,7 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
// that on the i-th row the diagonal is negative,
// multiply U from the left by an identity matrix that
// has -1 on the i-th row.
blas64.Scal(n-i, -1, blas64.Vector{Data: umat.Data[i*stride+i : i*stride+n], Inc: 1})
blas64.Scal(-1, blas64.Vector{N: n - i, Data: umat.Data[i*stride+i : i*stride+n], Inc: 1})
}
}
if ok {
Expand Down

0 comments on commit 30781f0

Please sign in to comment.