Merge ffc4d79 into 39b486f

mat/band.go

@@ -0,0 +1,163 @@
+// 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 {
+	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) == r*(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.
+//
+// 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) != r*bc {
+		panic(ErrShape)
+	}
+	if data == nil {
+		data = make([]float64, r*bc)
+	}
+	return &BandDense{
+		mat: blas64.Band{
+			Rows:   r,
+			Cols:   c,
+			KL:     kl,
+			KU:     ku,
+			Stride: bc,
+			Data:   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
+}

mat/band_test.go

@@ -0,0 +1,207 @@
+// 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 (
+	"reflect"
+	"testing"
+
+	"gonum.org/v1/gonum/blas/blas64"
+)
+
+func TestNewBand(t *testing.T) {
+	for i, test := range []struct {
+		data   []float64
+		r, c   int
+		kl, ku int
+		mat    *BandDense
+		dense  *Dense
+	}{
+		{
+			data: []float64{
+				-1, 1, 2, 3,
+				4, 5, 6, 7,
+				8, 9, 10, 11,
+				12, 13, 14, 15,
+				16, 17, 18, -1,
+				19, 20, -1, -1,
+			},
+			r: 6, c: 6,
+			kl: 1, ku: 2,
+			mat: &BandDense{
+				mat: blas64.Band{
+					Rows:   6,
+					Cols:   6,
+					KL:     1,
+					KU:     2,
+					Stride: 4,
+					Data: []float64{
+						-1, 1, 2, 3,
+						4, 5, 6, 7,
+						8, 9, 10, 11,
+						12, 13, 14, 15,
+						16, 17, 18, -1,
+						19, 20, -1, -1,
+					},
+				},
+			},
+			dense: NewDense(6, 6, []float64{
+				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,
+			}),
+		},
+	} {
+		band := NewBandDense(test.r, test.c, test.kl, test.ku, test.data)
+		rows, cols := band.Dims()
+
+		if rows != test.r {
+			t.Errorf("unexpected number of rows for test %d: got: %d want: %d", i, rows, test.r)
+		}
+		if cols != test.c {
+			t.Errorf("unexpected number of cols for test %d: got: %d want: %d", i, cols, test.c)
+		}
+		if !reflect.DeepEqual(band, test.mat) {
+			t.Errorf("unexpected value via reflect for test %d: got: %v want: %v", i, band, test.mat)
+		}
+		if !Equal(band, test.mat) {
+			t.Errorf("unexpected value via mat.Equal for test %d: got: %v want: %v", i, band, test.mat)
+		}
+		if !Equal(band, test.dense) {
+			t.Errorf("unexpected value via mat.Equal(band, dense) for test %d:\ngot:\n% v\nwant:\n% v", i, Formatted(band), Formatted(test.dense))
+		}
+	}
+}
+
+func TestBandAtSet(t *testing.T) {
+	// 2  3  4  0  0  0
+	// 5  6  7  8  0  0
+	// 0  9 10 11 12  0
+	// 0  0 13 14 15 16
+	// 0  0  0 17 18 19
+	// 0  0  0  0 21 22
+	band := NewBandDense(6, 6, 1, 2, []float64{
+		-1, 2, 3, 4,
+		5, 6, 7, 8,
+		9, 10, 11, 12,
+		13, 14, 15, 16,
+		17, 18, 19, -1,
+		21, 22, -1, -1,
+	})
+
+	rows, cols := band.Dims()
+	kl, ku := band.Bandwidth()
+
+	// Explicitly test all indexes.
+	want := bandImplicit{rows, cols, kl, ku, func(i, j int) float64 {
+		return float64(i*(kl+ku) + j + kl + 1)
+	}}
+	for i := 0; i < 6; i++ {
+		for j := 0; j < 6; j++ {
+			if band.At(i, j) != want.At(i, j) {
+				t.Errorf("unexpected value for band.At(%d, %d): got:%v want:%v", band.At(i, j), want.At(i, j))
+			}
+		}
+	}
+	// Do that same thing via a call to Equal.
+	if !Equal(band, want) {
+		t.Errorf("unexpected value via mat.Equal: got: %v want: %v", band, want)
+	}
+
+	// Check At out of bounds
+	for _, row := range []int{-1, rows, rows + 1} {
+		panicked, message := panics(func() { band.At(row, 0) })
+		if !panicked || message != ErrRowAccess.Error() {
+			t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
+		}
+	}
+	for _, col := range []int{-1, cols, cols + 1} {
+		panicked, message := panics(func() { band.At(0, col) })
+		if !panicked || message != ErrColAccess.Error() {
+			t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
+		}
+	}
+
+	// Check Set out of bounds
+	for _, row := range []int{-1, rows, rows + 1} {
+		panicked, message := panics(func() { band.SetBand(row, 0, 1.2) })
+		if !panicked || message != ErrRowAccess.Error() {
+			t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
+		}
+	}
+	for _, col := range []int{-1, cols, cols + 1} {
+		panicked, message := panics(func() { band.SetBand(0, col, 1.2) })
+		if !panicked || message != ErrColAccess.Error() {
+			t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
+		}
+	}
+
+	for _, st := range []struct {
+		row, col int
+	}{
+		{row: 0, col: 3},
+		{row: 0, col: 4},
+		{row: 0, col: 5},
+		{row: 1, col: 4},
+		{row: 1, col: 5},
+		{row: 2, col: 5},
+		{row: 2, col: 0},
+		{row: 3, col: 1},
+		{row: 4, col: 2},
+		{row: 5, col: 3},
+	} {
+		panicked, message := panics(func() { band.SetBand(st.row, st.col, 1.2) })
+		if !panicked || message != ErrBandSet.Error() {
+			t.Errorf("expected panic for %+v %s", st, message)
+		}
+	}
+
+	for _, st := range []struct {
+		row, col  int
+		orig, new float64
+	}{
+		{row: 1, col: 2, orig: 7, new: 15},
+		{row: 2, col: 3, orig: 11, new: 15},
+	} {
+		if e := band.At(st.row, st.col); e != st.orig {
+			t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.row, st.col, e, st.orig)
+		}
+		band.SetBand(st.row, st.col, st.new)
+		if e := band.At(st.row, st.col); e != st.new {
+			t.Errorf("unexpected value for At(%d, %d) after SetTri(%[1]d, %d, %v): got: %v want: %[3]v", st.row, st.col, st.new, e)
+		}
+	}
+}
+
+// bandImplicit is an implicit band matrix returning val(i, j)
+// for the value at (i, j).
+type bandImplicit struct {
+	r, c, kl, ku int
+	val          func(i, j int) float64
+}
+
+func (b bandImplicit) Dims() (r, c int) {
+	return b.r, b.c
+}
+
+func (b bandImplicit) T() Matrix {
+	return Transpose{b}
+}
+
+func (b bandImplicit) At(i, j int) float64 {
+	if i < 0 || b.r <= i {
+		panic("row")
+	}
+	if j < 0 || b.c <= j {
+		panic("col")
+	}
+	if j < i-b.kl || i+b.ku < j {
+		return 0
+	}
+	return b.val(i, j)
+}

mat/errors.go

@@ -129,6 +129,7 @@ var (
 	ErrPivot               = Error{"matrix: malformed pivot list"}
 	ErrTriangle            = Error{"matrix: triangular storage mismatch"}
 	ErrTriangleSet         = Error{"matrix: triangular set out of bounds"}
+	ErrBandSet             = Error{"matrix: band set out of bounds"}
 	ErrSliceLengthMismatch = Error{"matrix: input slice length mismatch"}
 	ErrNotPSD              = Error{"matrix: input not positive symmetric definite"}
 	ErrFailedEigen         = Error{"matrix: eigendecomposition not successful"}