Skip to content
This repository has been archived by the owner on Dec 10, 2018. It is now read-only.

Commit

Permalink
Change Triangular to TriDense and add Triangular interface.
Browse files Browse the repository at this point in the history 
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
  • Loading branch information
@btracey
btracey committed Mar 26, 2015
1 parent 2bf9d4e commit 8a38c9f
Show file tree
Hide file tree
Showing 6 changed files with 102 additions and 44 deletions.
30 changes: 19 additions & 11 deletions mat64/cholesky.go
View file
Expand Up @@ -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{
Expand Down Expand Up @@ -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 {
Expand All @@ -93,22 +93,26 @@ 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)
}
}

// 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)
Expand All @@ -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)
}
Expand Down
22 changes: 11 additions & 11 deletions mat64/cholesky_test.go
View file
Expand Up @@ -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
}{
{
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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)

Expand Down
8 changes: 4 additions & 4 deletions mat64/index_bound_checks.go
View file
Expand Up @@ -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)
}
Expand All @@ -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)
}
Expand Down
8 changes: 4 additions & 4 deletions mat64/index_no_bound_checks.go
View file
Expand Up @@ -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)
}
Expand All @@ -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
Expand All @@ -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)
}
Expand All @@ -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
}
70 changes: 60 additions & 10 deletions mat64/triangular.go
View file
Expand Up @@ -6,23 +6,37 @@ 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
// mat == nil, new data will be allocated, and will panic if neither of these
// 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")
}
Expand All @@ -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,
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
}
8 changes: 4 additions & 4 deletions mat64/triangular_test.go
View file
Expand Up @@ -11,7 +11,7 @@ func (s *S) TestNewTriangular(c *check.C) {
data []float64
N int
upper bool
mat *Triangular
mat *TriDense
}{
{
data: []float64{
Expand All @@ -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,
Expand All @@ -30,15 +30,15 @@ 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))
c.Check(t, check.DeepEquals, test.mat, check.Commentf("Test %d", i))
}
}
func (s *S) TestTriAtSet(c *check.C) {
t := &Triangular{blas64.Triangular{
t := &TriDense{blas64.Triangular{
N: 3,
Stride: 3,
Uplo: blas.Upper,
Expand Down

0 comments on commit 8a38c9f

Please sign in to comment.