From 8a38c9f4a32319708676835f535719faf6919360 Mon Sep 17 00:00:00 2001 From: btracey Date: Thu, 26 Mar 2015 00:11:14 -0700 Subject: [PATCH] Change Triangular to TriDense and add Triangular interface. This brings the Triangular type in line with Dense and Symmetric Made SolveChol take in a Triangular shortened Triangular path Changed method name to Triangle --- mat64/cholesky.go | 30 +++++++++------ mat64/cholesky_test.go | 22 +++++------ mat64/index_bound_checks.go | 8 ++-- mat64/index_no_bound_checks.go | 8 ++-- mat64/triangular.go | 70 +++++++++++++++++++++++++++++----- mat64/triangular_test.go | 8 ++-- 6 files changed, 102 insertions(+), 44 deletions(-) diff --git a/mat64/cholesky.go b/mat64/cholesky.go index e5ac5ea..7e5a2be 100644 --- a/mat64/cholesky.go +++ b/mat64/cholesky.go @@ -19,7 +19,7 @@ const badTriangle = "mat64: invalid triangle" // whether the matrix is positive definite. The returned matrix is either a // lower triangular matrix such that A = L * L^T or an upper triangular matrix // such that A = U^T * U depending on the upper parameter. -func (t *Triangular) Cholesky(a *SymDense, upper bool) (ok bool) { +func (t *TriDense) Cholesky(a *SymDense, upper bool) (ok bool) { n := a.Symmetric() if t.isZero() { t.mat = blas64.Triangular{ @@ -81,7 +81,7 @@ func (t *Triangular) Cholesky(a *SymDense, upper bool) (ok bool) { // SolveCholesky finds the matrix x that solves A * X = B where A = L * L^T or // A = U^T * U, and U or L are represented by t. The matrix A must be symmetric // and positive definite. -func (m *Dense) SolveCholesky(t *Triangular, b Matrix) { +func (m *Dense) SolveCholesky(t Triangular, b Matrix) { _, n := t.Dims() bm, bn := b.Dims() if n != bm { @@ -93,13 +93,17 @@ func (m *Dense) SolveCholesky(t *Triangular, b Matrix) { m.Copy(b) } - switch t.mat.Uplo { + // TODO(btracey): Implement an algorithm that doesn't require a copy into + // a blas64.Triangular. + ta := getBlasTriangular(t) + + switch ta.Uplo { case blas.Upper: - blas64.Trsm(blas.Left, blas.Trans, 1, t.mat, m.mat) - blas64.Trsm(blas.Left, blas.NoTrans, 1, t.mat, m.mat) + blas64.Trsm(blas.Left, blas.Trans, 1, ta, m.mat) + blas64.Trsm(blas.Left, blas.NoTrans, 1, ta, m.mat) case blas.Lower: - blas64.Trsm(blas.Left, blas.NoTrans, 1, t.mat, m.mat) - blas64.Trsm(blas.Left, blas.Trans, 1, t.mat, m.mat) + blas64.Trsm(blas.Left, blas.NoTrans, 1, ta, m.mat) + blas64.Trsm(blas.Left, blas.Trans, 1, ta, m.mat) default: panic(badTriangle) } @@ -107,8 +111,8 @@ func (m *Dense) SolveCholesky(t *Triangular, b Matrix) { // SolveTri finds the matrix x that solves op(A) * X = B where A is a triangular // matrix and op is specified by trans. -func (m *Dense) SolveTri(a *Triangular, trans bool, b Matrix) { - _, n := a.Dims() +func (m *Dense) SolveTri(a Triangular, trans bool, b Matrix) { + n, _ := a.Triangle() bm, bn := b.Dims() if n != bm { panic(ErrShape) @@ -119,13 +123,17 @@ func (m *Dense) SolveTri(a *Triangular, trans bool, b Matrix) { m.Copy(b) } + // TODO(btracey): Implement an algorithm that doesn't require a copy into + // a blas64.Triangular. + ta := getBlasTriangular(a) + t := blas.NoTrans if trans { t = blas.Trans } - switch a.mat.Uplo { + switch ta.Uplo { case blas.Upper, blas.Lower: - blas64.Trsm(blas.Left, t, 1, a.mat, m.mat) + blas64.Trsm(blas.Left, t, 1, ta, m.mat) default: panic(badTriangle) } diff --git a/mat64/cholesky_test.go b/mat64/cholesky_test.go index 45d963c..386da8e 100644 --- a/mat64/cholesky_test.go +++ b/mat64/cholesky_test.go @@ -10,9 +10,9 @@ func (s *S) TestCholesky(c *check.C) { for _, t := range []struct { a *SymDense upper bool - f *Triangular + f *TriDense - want *Triangular + want *TriDense pd bool }{ { @@ -22,9 +22,9 @@ func (s *S) TestCholesky(c *check.C) { 0, 0, 6, }), upper: false, - f: &Triangular{}, + f: &TriDense{}, - want: NewTriangular(3, false, []float64{ + want: NewTriDense(3, false, []float64{ 2, 0, 0, 0.5, 1.3228756555322954, 0, 0.5, 2.0788046015507495, 1.195228609334394, @@ -38,9 +38,9 @@ func (s *S) TestCholesky(c *check.C) { 0, 0, 6, }), upper: true, - f: &Triangular{}, + f: &TriDense{}, - want: NewTriangular(3, true, []float64{ + want: NewTriDense(3, true, []float64{ 2, 0.5, 0.5, 0, 1.3228756555322954, 2.0788046015507495, 0, 0, 1.195228609334394, @@ -54,9 +54,9 @@ func (s *S) TestCholesky(c *check.C) { 0, 0, 6, }), upper: false, - f: NewTriangular(3, false, nil), + f: NewTriDense(3, false, nil), - want: NewTriangular(3, false, []float64{ + want: NewTriDense(3, false, []float64{ 2, 0, 0, 0.5, 1.3228756555322954, 0, 0.5, 2.0788046015507495, 1.195228609334394, @@ -70,9 +70,9 @@ func (s *S) TestCholesky(c *check.C) { 0, 0, 6, }), upper: true, - f: NewTriangular(3, false, nil), + f: NewTriDense(3, false, nil), - want: NewTriangular(3, true, []float64{ + want: NewTriDense(3, true, []float64{ 2, 0.5, 0.5, 0, 1.3228756555322954, 2.0788046015507495, 0, 0, 1.195228609334394, @@ -126,7 +126,7 @@ func (s *S) TestCholeskySolve(c *check.C) { ans: NewDense(2, 1, []float64{5, 6}), }, } { - var f Triangular + var f TriDense ok := f.Cholesky(t.a, false) c.Assert(ok, check.Equals, true) diff --git a/mat64/index_bound_checks.go b/mat64/index_bound_checks.go index e4a5e18..5de12af 100644 --- a/mat64/index_bound_checks.go +++ b/mat64/index_bound_checks.go @@ -103,11 +103,11 @@ func (t *SymDense) set(r, c int, v float64) { } // At returns the element at row r and column c. -func (t *Triangular) At(r, c int) float64 { +func (t *TriDense) At(r, c int) float64 { return t.at(r, c) } -func (t *Triangular) at(r, c int) float64 { +func (t *TriDense) at(r, c int) float64 { if r >= t.mat.N || r < 0 { panic(ErrRowAccess) } @@ -128,11 +128,11 @@ func (t *Triangular) at(r, c int) float64 { // SetTri sets the element of the triangular matrix at row r and column c. // Set panics if the location is outside the appropriate half of the matrix. -func (t *Triangular) SetTri(r, c int, v float64) { +func (t *TriDense) SetTri(r, c int, v float64) { t.set(r, c, v) } -func (t *Triangular) set(r, c int, v float64) { +func (t *TriDense) set(r, c int, v float64) { if r >= t.mat.N || r < 0 { panic(ErrRowAccess) } diff --git a/mat64/index_no_bound_checks.go b/mat64/index_no_bound_checks.go index f9b89a6..27f880e 100644 --- a/mat64/index_no_bound_checks.go +++ b/mat64/index_no_bound_checks.go @@ -103,7 +103,7 @@ func (s *SymDense) set(r, c int, v float64) { } // At returns the element at row r and column c. -func (t *Triangular) At(r, c int) float64 { +func (t *TriDense) At(r, c int) float64 { if r >= t.mat.N || r < 0 { panic(ErrRowAccess) } @@ -113,7 +113,7 @@ func (t *Triangular) At(r, c int) float64 { return t.at(r, c) } -func (t *Triangular) at(r, c int) float64 { +func (t *TriDense) at(r, c int) float64 { if t.mat.Uplo == blas.Upper { if r > c { return 0 @@ -128,7 +128,7 @@ func (t *Triangular) at(r, c int) float64 { // SetTri sets the element at row r and column c. Set panics if the location is outside // the appropriate half of the matrix. -func (t *Triangular) SetTri(r, c int, v float64) { +func (t *TriDense) SetTri(r, c int, v float64) { if r >= t.mat.N || r < 0 { panic(ErrRowAccess) } @@ -144,6 +144,6 @@ func (t *Triangular) SetTri(r, c int, v float64) { t.set(r, c, v) } -func (t *Triangular) set(r, c int, v float64) { +func (t *TriDense) set(r, c int, v float64) { t.mat.Data[r*t.mat.Stride+c] = v } diff --git a/mat64/triangular.go b/mat64/triangular.go index 9286f83..3364a51 100644 --- a/mat64/triangular.go +++ b/mat64/triangular.go @@ -6,15 +6,29 @@ import ( ) var ( - triangular *Triangular - _ Matrix = triangular + triDense *TriDense + _ Matrix = triDense + _ Triangular = triDense + _ RawTriangular = triDense ) -// Triangular represents an upper or lower triangular matrix. -type Triangular struct { +// TriDense represents an upper or lower triangular matrix in dense storage +// format. +type TriDense struct { mat blas64.Triangular } +type Triangular interface { + Matrix + // Triangular returns the number of rows/columns in the matrix and if it is + // an upper triangular matrix. + Triangle() (n int, upper bool) +} + +type RawTriangular interface { + RawTriangular() blas64.Triangular +} + // NewTriangular constructs an n x n triangular matrix. The constructed matrix // is upper triangular if upper == true and lower triangular otherwise. // If len(mat) == n * n, mat will be used to hold the underlying data, if @@ -22,7 +36,7 @@ type Triangular struct { // cases is true. // The underlying data representation is the same as that of a Dense matrix, // except the values of the entries in the opposite half are completely ignored. -func NewTriangular(n int, upper bool, mat []float64) *Triangular { +func NewTriDense(n int, upper bool, mat []float64) *TriDense { if n < 0 { panic("mat64: negative dimension") } @@ -36,7 +50,7 @@ func NewTriangular(n int, upper bool, mat []float64) *Triangular { if upper { uplo = blas.Upper } - return &Triangular{blas64.Triangular{ + return &TriDense{blas64.Triangular{ N: n, Stride: n, Data: mat, @@ -45,15 +59,19 @@ func NewTriangular(n int, upper bool, mat []float64) *Triangular { }} } -func (t *Triangular) Dims() (r, c int) { +func (t *TriDense) Dims() (r, c int) { return t.mat.N, t.mat.N } -func (t *Triangular) RawTriangular() blas64.Triangular { +func (t *TriDense) Triangle() (n int, upper bool) { + return t.mat.N, t.mat.Uplo == blas.Upper +} + +func (t *TriDense) RawTriangular() blas64.Triangular { return t.mat } -func (t *Triangular) isZero() bool { +func (t *TriDense) isZero() bool { // It must be the case that t.Dims() returns // zeros in this case. See comment in Reset(). return t.mat.Stride == 0 @@ -63,7 +81,7 @@ func (t *Triangular) isZero() bool { // receiver of a dimensionally restricted operation. // // See the Reseter interface for more information. -func (t *Triangular) Reset() { +func (t *TriDense) Reset() { // No change of Stride, N to 0 may // be made unless both are set to 0. t.mat.N, t.mat.Stride = 0, 0 @@ -72,3 +90,35 @@ func (t *Triangular) Reset() { t.mat.Uplo = 0 t.mat.Data = t.mat.Data[:0] } + +// getBlasTriangular transforms t into a blas64.Triangular. If t is a RawTriangular, +// the direct matrix representation is returned, otherwise t is copied into one. +func getBlasTriangular(t Triangular) blas64.Triangular { + n, upper := t.Triangle() + rt, ok := t.(RawTriangular) + if ok { + return rt.RawTriangular() + } + ta := blas64.Triangular{ + N: n, + Stride: n, + Diag: blas.NonUnit, + Data: make([]float64, n*n), + } + if upper { + ta.Uplo = blas.Upper + for i := 0; i < n; i++ { + for j := i; j < n; j++ { + ta.Data[i*n+j] = t.At(i, j) + } + } + return ta + } + ta.Uplo = blas.Lower + for i := 0; i < n; i++ { + for j := 0; j < i; j++ { + ta.Data[i*n+j] = t.At(i, j) + } + } + return ta +} diff --git a/mat64/triangular_test.go b/mat64/triangular_test.go index 7093012..786819d 100644 --- a/mat64/triangular_test.go +++ b/mat64/triangular_test.go @@ -11,7 +11,7 @@ func (s *S) TestNewTriangular(c *check.C) { data []float64 N int upper bool - mat *Triangular + mat *TriDense }{ { data: []float64{ @@ -21,7 +21,7 @@ func (s *S) TestNewTriangular(c *check.C) { }, N: 3, upper: true, - mat: &Triangular{blas64.Triangular{ + mat: &TriDense{blas64.Triangular{ N: 3, Stride: 3, Uplo: blas.Upper, @@ -30,7 +30,7 @@ func (s *S) TestNewTriangular(c *check.C) { }}, }, } { - t := NewTriangular(test.N, test.upper, test.data) + t := NewTriDense(test.N, test.upper, test.data) rows, cols := t.Dims() c.Check(rows, check.Equals, test.N, check.Commentf("Test %d", i)) c.Check(cols, check.Equals, test.N, check.Commentf("Test %d", i)) @@ -38,7 +38,7 @@ func (s *S) TestNewTriangular(c *check.C) { } } func (s *S) TestTriAtSet(c *check.C) { - t := &Triangular{blas64.Triangular{ + t := &TriDense{blas64.Triangular{ N: 3, Stride: 3, Uplo: blas.Upper,