Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

mat: add basis for band matrix type #95

Merged
merged 3 commits into from Jul 8, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
170 changes: 170 additions & 0 deletions mat/band.go
@@ -0,0 +1,170 @@
// Copyright ©2017 The gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package mat

import (
"gonum.org/v1/gonum/blas/blas64"
)

var (
bandDense *BandDense
_ Matrix = bandDense
_ Banded = bandDense
_ RawBand = bandDense
)

// BandDense represents a band matrix in dense storage format.
type BandDense struct {
mat blas64.Band
}

type Banded interface {
Matrix
// Bandwidth returns the lower and upper bandwidth values for
// the matrix. The total bandwidth of the matrix is kl+ku+1.
Bandwidth() (kl, ku int)

// TBand is the equivalent of the T() method in the Matrix
// interface but guarantees the transpose is of banded type.
TBand() Banded
}

type RawBand interface {
RawBand() blas64.Band
}

var (
_ Matrix = TransposeBand{}
_ Banded = TransposeBand{}
_ UntransposeBander = TransposeBand{}
)

// TransposeBand is a type for performing an implicit transpose of a band
// matrix. It implements the Banded interface, returning values from the
// transpose of the matrix within.
type TransposeBand struct {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be TransposeBanded? In general, from a quick look, I don't understand the rules when Band or Banded is used.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did write it like that, but reverted because of length. Various people call these matrices band or banded. Here Banded is the interface and Band is the theoretical notion. I don't really mind either way.

Banded Banded
}

// At returns the value of the element at row i and column j of the transposed
// matrix, that is, row j and column i of the Banded field.
func (t TransposeBand) At(i, j int) float64 {
return t.Banded.At(j, i)
}

// Dims returns the dimensions of the transposed matrix.
func (t TransposeBand) Dims() (r, c int) {
c, r = t.Banded.Dims()
return r, c
}

// T performs an implicit transpose by returning the Banded field.
func (t TransposeBand) T() Matrix {
return t.Banded
}

// Bandwidth returns the number of rows/columns in the matrix and its orientation.
func (t TransposeBand) Bandwidth() (kl, ku int) {
kl, ku = t.Banded.Bandwidth()
return ku, kl
}

// TBand performs an implicit transpose by returning the Banded field.
func (t TransposeBand) TBand() Banded {
return t.Banded
}

// Untranspose returns the Banded field.
func (t TransposeBand) Untranspose() Matrix {
return t.Banded
}

// UntransposeBand returns the Banded field.
func (t TransposeBand) UntransposeBand() Banded {
return t.Banded
}

// NewBandDense creates a new Band matrix with r rows and c columns. If data == nil,
// a new slice is allocated for the backing slice. If len(data) == min(r, c+kl)*(kl+ku+1),
// data is used as the backing slice, and changes to the elements of the returned
// BandDense will be reflected in data. If neither of these is true, NewBandDense
// will panic. kl must be at least zero and less r, and ku must be at least zero and
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note comments on kl and ku now.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this behavior match Lapack? I think we want to be compliant here. I've only started looking at lapack banded implementations, but it would be nice to not have to play too many games (that said, we are stricter for Vector and it's worked out fine)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The documentation says that the * elements are not referenced. If this causes a failure, I think then that the LAPACK implementation is buggy.

Looking at how we do this, I suspect it will fail for the gonum implementation, for example here. However, I think that is a result of our slice optimisations; the Fortran should not fail.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. I agree. We/I should probably revisit the banded BLAS implementations before we're 6 layers deep in Lapack trying to figure out what's going on.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are two ways to deal with those, either continue if the proposed length is 0, or probably preferably, make the same assessment that I am doing here.

// less than c, otherwise NewBandDense will panic.
//
// The data must be arranged in row-major order constructed by removing the zeros
// from the rows outside the band and aligning the diagonals. For example, the matrix
// 1 2 3 0 0 0
// 4 5 6 7 0 0
// 0 8 9 10 11 0
// 0 0 12 13 14 15
// 0 0 0 16 17 18
// 0 0 0 0 19 20
// becomes (* entries are never accessed)
// * 1 2 3
// 4 5 6 7
// 8 9 10 11
// 12 13 14 15
// 16 17 18 *
// 19 20 * *
// which is passed to NewBandDense as []float64{*, 1, 2, 3, 4, ...} with kl=1 and ku=2.
// Only the values in the band portion of the matrix are used.
func NewBandDense(r, c, kl, ku int, data []float64) *BandDense {
if r < 0 || c < 0 || kl < 0 || ku < 0 {
panic("mat: negative dimension")
}
if kl+1 > r || ku+1 > c {
panic("mat: band out of range")
}
bc := kl + ku + 1
if data != nil && len(data) != min(r, c+kl)*bc {
panic(ErrShape)
}
if data == nil {
data = make([]float64, min(r, c+kl)*bc)
}
return &BandDense{
mat: blas64.Band{
Rows: r,
Cols: c,
KL: kl,
KU: ku,
Stride: bc,
Data: data,
},
}
}

// NewDiagonalRect is a convenience function that returns a diagonal matrix represented by a
// BandDense. The length of data must be min(r, c) otherwise NewDiagonalRect will panic.
func NewDiagonalRect(r, c int, data []float64) *BandDense {
return NewBandDense(r, c, 0, 0, data)
}

// Dims returns the number of rows and columns in the matrix.
func (b *BandDense) Dims() (r, c int) {
return b.mat.Rows, b.mat.Cols
}

// Bandwidth returns the upper and lower bandwidths of the matrix.
func (b *BandDense) Bandwidth() (kl, ku int) {
return b.mat.KL, b.mat.KU
}

// T performs an implicit transpose by returning the receiver inside a Transpose.
func (b *BandDense) T() Matrix {
return Transpose{b}
}

// TBand performs an implicit transpose by returning the receiver inside a TransposeBand.
func (b *BandDense) TBand() Banded {
return TransposeBand{b}
}

// RawBand returns the underlying blas64.Band used by the receiver.
// Changes to elements in the receiver following the call will be reflected
// in returned blas64.Band.
func (b *BandDense) RawBand() blas64.Band {
return b.mat
}