blas/blas64: fix failures and reuse band logic
MODIFY
blas/blas64/conv.go
blas/blas64/conv.go
@@ -106,14 +106,18 @@ func (t Symmetric) From(a SymmetricCols) {
type TriangularCols Triangular
// From fills the receiver with elements from a. The receiver
-// must have the same dimensions as a and have adequate backing
-// data storage.
+// must have the same dimensions, uplo and diag as a and have
+// adequate backing data storage.
func (t TriangularCols) From(a Triangular) {
if t.N != a.N {
panic("blas64: mismatched dimension")
}
- t.Uplo = a.Uplo
- t.Diag = a.Diag
+ if t.Uplo != a.Uplo {
+ panic("blas64: mismatched BLAS uplo")
+ }
+ if t.Diag != a.Diag {
+ panic("blas64: mismatched BLAS diag")
+ }
switch a.Uplo {
default:
panic("blas64: bad BLAS uplo")
@@ -139,14 +143,18 @@ func (t TriangularCols) From(a Triangular) {
}
// From fills the receiver with elements from a. The receiver
-// must have the same dimensions as a and have adequate backing
-// data storage.
+// must have the same dimensions, uplo and diag as a and have
+// adequate backing data storage.
func (t Triangular) From(a TriangularCols) {
if t.N != a.N {
panic("blas64: mismatched dimension")
}
- t.Uplo = a.Uplo
- t.Diag = a.Diag
+ if t.Uplo != a.Uplo {
+ panic("blas64: mismatched BLAS uplo")
+ }
+ if t.Diag != a.Diag {
+ panic("blas64: mismatched BLAS diag")
+ }
switch a.Uplo {
default:
panic("blas64: bad BLAS uplo")
@@ -242,28 +250,78 @@ func (t SymmetricBandCols) From(a SymmetricBand) {
if t.Uplo != a.Uplo {
panic("blas64: mismatched BLAS uplo")
}
+ dst := BandCols{
+ Rows: t.N, Cols: t.N,
+ Stride: t.Stride,
+ Data: t.Data,
+ }
+ src := Band{
+ Rows: a.N, Cols: a.N,
+ Stride: a.Stride,
+ Data: a.Data,
+ }
switch a.Uplo {
default:
panic("blas64: bad BLAS uplo")
case blas.Upper:
- for i := 0; i < a.N; i++ {
- for j := i; j < min(i+a.K+1, a.N); j++ {
- t.Data[i+t.K-j+j*t.Stride] = a.Data[j+a.K-i+i*a.Stride]
+ dst.KU = t.K
+ src.KU = a.K
+ case blas.Lower:
+ dst.KL = t.K
+ src.KL = a.K
}
+ dst.From(src)
}
- case blas.Lower:
- for i := 0; i < a.N; i++ {
- for j := max(0, i-a.K); j <= i; j++ {
- t.Data[i+t.K-j+j*t.Stride] = a.Data[j+a.K-i+i*a.Stride]
+
+// From fills the receiver with elements from a. The receiver
+// must have the same dimensions, bandwidth and uplo as a and
+// have adequate backing data storage.
+func (t SymmetricBand) From(a SymmetricBandCols) {
+ if t.N != a.N {
+ panic("blas64: mismatched dimension")
+ }
+ if t.K != a.K {
+ panic("blas64: mismatched bandwidth")
}
+ if a.Stride < a.K+1 {
+ panic("blas64: short stride for source")
+ }
+ if t.Stride < t.K+1 {
+ panic("blas64: short stride for destination")
}
+ if t.Uplo != a.Uplo {
+ panic("blas64: mismatched BLAS uplo")
+ }
+ dst := Band{
+ Rows: t.N, Cols: t.N,
+ Stride: t.Stride,
+ Data: t.Data,
}
+ src := BandCols{
+ Rows: a.N, Cols: a.N,
+ Stride: a.Stride,
+ Data: a.Data,
+ }
+ switch a.Uplo {
+ default:
+ panic("blas64: bad BLAS uplo")
+ case blas.Upper:
+ dst.KU = t.K
+ src.KU = a.K
+ case blas.Lower:
+ dst.KL = t.K
+ src.KL = a.K
+ }
+ dst.From(src)
}
+// TriangularBandCols represents a symmetric matrix using the band column-major storage scheme.
+type TriangularBandCols TriangularBand
+
// From fills the receiver with elements from a. The receiver
// must have the same dimensions, bandwidth and uplo as a and
// have adequate backing data storage.
-func (t SymmetricBand) From(a SymmetricBandCols) {
+func (t TriangularBandCols) From(a TriangularBand) {
if t.N != a.N {
panic("blas64: mismatched dimension")
}
@@ -279,22 +337,75 @@ func (t SymmetricBand) From(a SymmetricBandCols) {
if t.Uplo != a.Uplo {
panic("blas64: mismatched BLAS uplo")
}
+ if t.Diag != a.Diag {
+ panic("blas64: mismatched BLAS diag")
+ }
+ dst := BandCols{
+ Rows: t.N, Cols: t.N,
+ Stride: t.Stride,
+ Data: t.Data,
+ }
+ src := Band{
+ Rows: a.N, Cols: a.N,
+ Stride: a.Stride,
+ Data: a.Data,
+ }
switch a.Uplo {
default:
panic("blas64: bad BLAS uplo")
case blas.Upper:
- for j := 0; j < a.N; j++ {
- for i := j; i < min(j+a.K+1, a.N); i++ {
- t.Data[j+a.K-i+i*a.Stride] = a.Data[i+t.K-j+j*t.Stride]
+ dst.KU = t.K
+ src.KU = a.K
+ case blas.Lower:
+ dst.KL = t.K
+ src.KL = a.K
}
+ dst.From(src)
}
- case blas.Lower:
- for j := 0; j < a.N; j++ {
- for i := max(0, j-a.K); i <= i; i++ {
- t.Data[j+a.K-i+i*a.Stride] = a.Data[i+t.K-j+j*t.Stride]
+
+// From fills the receiver with elements from a. The receiver
+// must have the same dimensions, bandwidth and uplo as a and
+// have adequate backing data storage.
+func (t TriangularBand) From(a TriangularBandCols) {
+ if t.N != a.N {
+ panic("blas64: mismatched dimension")
}
+ if t.K != a.K {
+ panic("blas64: mismatched bandwidth")
+ }
+ if a.Stride < a.K+1 {
+ panic("blas64: short stride for source")
}
+ if t.Stride < t.K+1 {
+ panic("blas64: short stride for destination")
+ }
+ if t.Uplo != a.Uplo {
+ panic("blas64: mismatched BLAS uplo")
+ }
+ if t.Diag != a.Diag {
+ panic("blas64: mismatched BLAS diag")
+ }
+ dst := Band{
+ Rows: t.N, Cols: t.N,
+ Stride: t.Stride,
+ Data: t.Data,
+ }
+ src := BandCols{
+ Rows: a.N, Cols: a.N,
+ Stride: a.Stride,
+ Data: a.Data,
+ }
+ switch a.Uplo {
+ default:
+ panic("blas64: bad BLAS uplo")
+ case blas.Upper:
+ dst.KU = t.K
+ src.KU = a.K
+ case blas.Lower:
+ dst.KL = t.K
+ src.KL = a.K
}
+ dst.From(src)
}
func min(a, b int) int {

MODIFY
blas/blas64/conv_test.go
blas/blas64/conv_test.go
@@ -499,7 +499,7 @@ var bandTests = []Band{
}
func TestConvertBand(t *testing.T) {
- for _, test := range bandTests[:1] {
+ for _, test := range bandTests {
colmajor := newBandColsFrom(test)
if !equalBand(colmajor, test) {
t.Errorf("unexpected result for row major to col major conversion:\n\tgot: %#v\n\tfrom:%#v",
@@ -530,14 +530,26 @@ func newSymmetricBandFrom(a SymmetricBandCols) SymmetricBand {
func (m SymmetricBand) n() (n int) { return m.N }
func (m SymmetricBand) at(i, j int) float64 {
+ b := Band{
+ Rows: m.N, Cols: m.N,
+ Stride: m.Stride,
+ Data: m.Data,
+ }
+ switch m.Uplo {
+ default:
+ panic("blas64: bad BLAS uplo")
+ case blas.Upper:
+ b.KU = m.K
if i > j {
i, j = j, i
}
- pj := j - i
- if m.K+1 <= pj {
- return 0
+ case blas.Lower:
+ b.KL = m.K
+ if i < j {
+ i, j = j, i
}
- return m.Data[i*m.Stride+pj]
+ }
+ return b.at(i, j)
}
func (m SymmetricBand) bandwidth() (k int) { return m.K }
func (m SymmetricBand) uplo() blas.Uplo { return m.Uplo }
@@ -559,14 +571,26 @@ func newSymmetricBandColsFrom(a SymmetricBand) SymmetricBandCols {
func (m SymmetricBandCols) n() (n int) { return m.N }
func (m SymmetricBandCols) at(i, j int) float64 {
+ b := BandCols{
+ Rows: m.N, Cols: m.N,
+ Stride: m.Stride,
+ Data: m.Data,
+ }
+ switch m.Uplo {
+ default:
+ panic("blas64: bad BLAS uplo")
+ case blas.Upper:
+ b.KU = m.K
if i > j {
i, j = j, i
}
- pj := j - i
- if m.K+1 <= pj {
- return 0
+ case blas.Lower:
+ b.KL = m.K
+ if i < j {
+ i, j = j, i
}
- return m.Data[j*m.Stride+pj]
+ }
+ return b.at(i, j)
}
func (m SymmetricBandCols) bandwidth() (k int) { return m.K }
func (m SymmetricBandCols) uplo() blas.Uplo { return m.Uplo }
@@ -667,7 +691,7 @@ var symmetricBandTests = []SymmetricBand{
}
func TestConvertSymBand(t *testing.T) {
- for _, test := range symmetricBandTests[:1] {
+ for _, test := range symmetricBandTests {
colmajor := newSymmetricBandColsFrom(test)
if !equalSymmetricBand(colmajor, test) {
t.Errorf("unexpected result for row major to col major conversion:\n\tgot: %#v\n\tfrom:%#v",
@@ -680,3 +704,209 @@ func TestConvertSymBand(t *testing.T) {
}
}
}
+
+func newTriangularBandFrom(a TriangularBandCols) TriangularBand {
+ t := TriangularBand{
+ N: a.N,
+ K: a.K,
+ Stride: a.K + 1,
+ Data: make([]float64, a.N*(a.K+1)),
+ Uplo: a.Uplo,
+ Diag: a.Diag,
+ }
+ for i := range t.Data {
+ t.Data[i] = math.NaN()
+ }
+ t.From(a)
+ return t
+}
+
+func (m TriangularBand) n() (n int) { return m.N }
+func (m TriangularBand) at(i, j int) float64 {
+ if m.Diag == blas.Unit && i == j {
+ return 1
+ }
+ b := Band{
+ Rows: m.N, Cols: m.N,
+ Stride: m.Stride,
+ Data: m.Data,
+ }
+ switch m.Uplo {
+ default:
+ panic("blas64: bad BLAS uplo")
+ case blas.Upper:
+ if i > j {
+ return 0
+ }
+ b.KU = m.K
+ case blas.Lower:
+ if i < j {
+ return 0
+ }
+ b.KL = m.K
+ }
+ return b.at(i, j)
+}
+func (m TriangularBand) bandwidth() (k int) { return m.K }
+func (m TriangularBand) uplo() blas.Uplo { return m.Uplo }
+func (m TriangularBand) diag() blas.Diag { return m.Diag }
+
+func newTriangularBandColsFrom(a TriangularBand) TriangularBandCols {
+ t := TriangularBandCols{
+ N: a.N,
+ K: a.K,
+ Stride: a.K + 1,
+ Data: make([]float64, a.N*(a.K+1)),
+ Uplo: a.Uplo,
+ Diag: a.Diag,
+ }
+ for i := range t.Data {
+ t.Data[i] = math.NaN()
+ }
+ t.From(a)
+ return t
+}
+
+func (m TriangularBandCols) n() (n int) { return m.N }
+func (m TriangularBandCols) at(i, j int) float64 {
+ if m.Diag == blas.Unit && i == j {
+ return 1
+ }
+ b := BandCols{
+ Rows: m.N, Cols: m.N,
+ Stride: m.Stride,
+ Data: m.Data,
+ }
+ switch m.Uplo {
+ default:
+ panic("blas64: bad BLAS uplo")
+ case blas.Upper:
+ if i > j {
+ return 0
+ }
+ b.KU = m.K
+ case blas.Lower:
+ if i < j {
+ return 0
+ }
+ b.KL = m.K
+ }
+ return b.at(i, j)
+}
+func (m TriangularBandCols) bandwidth() (k int) { return m.K }
+func (m TriangularBandCols) uplo() blas.Uplo { return m.Uplo }
+func (m TriangularBandCols) diag() blas.Diag { return m.Diag }
+
+type triangularBand interface {
+ n() (n int)
+ at(i, j int) float64
+ bandwidth() (k int)
+ uplo() blas.Uplo
+ diag() blas.Diag
+}
+
+func equalTriangularBand(a, b triangularBand) bool {
+ an := a.n()
+ bn := b.n()
+ if an != bn {
+ return false
+ }
+ if a.uplo() != b.uplo() {
+ return false
+ }
+ if a.diag() != b.diag() {
+ return false
+ }
+ ak := a.bandwidth()
+ bk := b.bandwidth()
+ if ak != bk {
+ return false
+ }
+ for i := 0; i < an; i++ {
+ for j := 0; j < an; j++ {
+ if a.at(i, j) != b.at(i, j) || math.IsNaN(a.at(i, j)) != math.IsNaN(b.at(i, j)) {
+ return false
+ }
+ }
+ }
+ return true
+}
+
+var triangularBandTests = []TriangularBand{
+ {N: 3, K: 0, Stride: 1, Uplo: blas.Upper, Data: []float64{
+ 1,
+ 2,
+ 3,
+ }},
+ {N: 3, K: 0, Stride: 1, Uplo: blas.Lower, Data: []float64{
+ 1,
+ 2,
+ 3,
+ }},
+ {N: 3, K: 1, Stride: 2, Uplo: blas.Upper, Data: []float64{
+ 1, 2,
+ 3, 4,
+ 5, -1,
+ }},
+ {N: 3, K: 1, Stride: 2, Uplo: blas.Lower, Data: []float64{
+ -1, 1,
+ 2, 3,
+ 4, 5,
+ }},
+ {N: 3, K: 2, Stride: 3, Uplo: blas.Upper, Data: []float64{
+ 1, 2, 3,
+ 4, 5, -1,
+ 6, -2, -3,
+ }},
+ {N: 3, K: 2, Stride: 3, Uplo: blas.Lower, Data: []float64{
+ -2, -1, 1,
+ -3, 2, 4,
+ 3, 5, 6,
+ }},
+
+ {N: 3, K: 0, Stride: 5, Uplo: blas.Upper, Data: []float64{
+ 1, 0, 0, 0, 0,
+ 2, 0, 0, 0, 0,
+ 3, 0, 0, 0, 0,
+ }},
+ {N: 3, K: 0, Stride: 5, Uplo: blas.Lower, Data: []float64{
+ 1, 0, 0, 0, 0,
+ 2, 0, 0, 0, 0,
+ 3, 0, 0, 0, 0,
+ }},
+ {N: 3, K: 1, Stride: 5, Uplo: blas.Upper, Data: []float64{
+ 1, 2, 0, 0, 0,
+ 3, 4, 0, 0, 0,
+ 5, -1, 0, 0, 0,
+ }},
+ {N: 3, K: 1, Stride: 5, Uplo: blas.Lower, Data: []float64{
+ -1, 1, 0, 0, 0,
+ 2, 3, 0, 0, 0,
+ 4, 5, 0, 0, 0,
+ }},
+ {N: 3, K: 2, Stride: 5, Uplo: blas.Upper, Data: []float64{
+ 1, 2, 3, 0, 0,
+ 4, 5, -1, 0, 0,
+ 6, -2, -3, 0, 0,
+ }},
+ {N: 3, K: 2, Stride: 5, Uplo: blas.Lower, Data: []float64{
+ -2, -1, 1, 0, 0,
+ -3, 2, 4, 0, 0,
+ 3, 5, 6, 0, 0,
+ }},
+}
+
+func TestConvertTriBand(t *testing.T) {
+ for _, test := range triangularBandTests {
+ colmajor := newTriangularBandColsFrom(test)
+ if !equalTriangularBand(colmajor, test) {
+ t.Errorf("unexpected result for row major to col major conversion:\n\tgot: %#v\n\tfrom:%#v",
+ colmajor, test)
+ }
+ rowmajor := newTriangularBandFrom(colmajor)
+ if !equalTriangularBand(rowmajor, test) {
+ t.Errorf("unexpected result for col major to row major conversion:\n\tgot: %#v\n\twant:%#v",
+ rowmajor, test)
+ }
+ }
+}

blas/blas64: fix failures and reuse band logic