diff --git a/referenceblas/level2double.go b/referenceblas/level2double.go index 85bbf462..bb9aac05 100644 --- a/referenceblas/level2double.go +++ b/referenceblas/level2double.go @@ -139,6 +139,83 @@ func (b Blas) Dgemv(o blas.Order, tA blas.Transpose, m, n int, alpha float64, a } } +// Dger performs the rank 1 operation +// A := alpha*x*y**T + A, +// where alpha is a scalar, x is an m element vector, y is an n element +// vector and A is an m by n matrix. +func (Blas) Dger(o blas.Order, m, n int, alpha float64, x []float64, incX int, y []float64, incY int, a []float64, lda int) { + // Check inputs + if o != blas.RowMajor && o != blas.ColMajor { + panic(badOrder) + } + if m < 0 { + panic("m < 0") + } + if n < 0 { + panic(negativeN) + } + if incX == 0 { + panic(zeroInc) + } + if incY == 0 { + panic(zeroInc) + } + if o == blas.RowMajor { + if lda > 1 && lda > n { + panic(badLda) + } + } else { + if lda > 1 && lda > m { + panic(badLda) + } + } + // Quick return if possible + if m == 0 || n == 0 || alpha == 0 { + return + } + + var jy, kx int + if incY > 0 { + jy = 1 + } else { + jy = -(n - 1) * incY + } + + if incY > 0 { + kx = 1 + } else { + kx = -(m - 1) * incX + } + + switch o { + default: + panic("should not be here") + case blas.RowMajor: + // TODO: Switch this to looping the other way + for j := 0; j < n; j++ { + if y[jy] != 0 { + tmp := alpha * y[jy] + ix := kx + for i := 0; i < m; i++ { + a[i*lda+j] += x[ix] * tmp + } + } + jy += incY + } + case blas.ColMajor: + for j := 0; j < n; j++ { + if y[jy] != 0 { + tmp := alpha * y[jy] + ix := kx + for i := 0; i < m; i++ { + a[j*lda+i] += x[ix] * tmp + } + } + jy += incY + } + } +} + // DTRMV performs one of the matrix-vector operations // x := A*x, or x := A**T*x, // where x is an n element vector and A is an n by n unit, or non-unit, @@ -448,83 +525,6 @@ func (b Blas) Dsymv(o blas.Order, ul blas.Uplo, n int, alpha float64, a []float6 } } -// Dger performs the rank 1 operation -// A := alpha*x*y**T + A, -// where alpha is a scalar, x is an m element vector, y is an n element -// vector and A is an m by n matrix. -func (Blas) Dger(o blas.Order, m, n int, alpha float64, x []float64, incX int, y []float64, incY int, a []float64, lda int) { - // Check inputs - if o != blas.RowMajor && o != blas.ColMajor { - panic(badOrder) - } - if m < 0 { - panic("m < 0") - } - if n < 0 { - panic(negativeN) - } - if incX == 0 { - panic(zeroInc) - } - if incY == 0 { - panic(zeroInc) - } - if o == blas.RowMajor { - if lda > 1 && lda > n { - panic(badLda) - } - } else { - if lda > 1 && lda > m { - panic(badLda) - } - } - // Quick return if possible - if m == 0 || n == 0 || alpha == 0 { - return - } - - var jy, kx int - if incY > 0 { - jy = 1 - } else { - jy = -(n - 1) * incY - } - - if incY > 0 { - kx = 1 - } else { - kx = -(m - 1) * incX - } - - switch o { - default: - panic("should not be here") - case blas.RowMajor: - // TODO: Switch this to looping the other way - for j := 0; j < n; j++ { - if y[jy] != 0 { - tmp := alpha * y[jy] - ix := kx - for i := 0; i < m; i++ { - a[i*lda+j] += x[ix] * tmp - } - } - jy += incY - } - case blas.ColMajor: - for j := 0; j < n; j++ { - if y[jy] != 0 { - tmp := alpha * y[jy] - ix := kx - for i := 0; i < m; i++ { - a[j*lda+i] += x[ix] * tmp - } - } - jy += incY - } - } -} - /* // Level 2 routines. Dgbmv(o Order, tA Transpose, m, n, kL, kU int, alpha float64, a []float64, lda int, x []float64, incX int, beta float64, y []float64, incY int) diff --git a/testblas/aux.go b/testblas/aux.go new file mode 100644 index 00000000..5b78d825 --- /dev/null +++ b/testblas/aux.go @@ -0,0 +1,79 @@ +package testblas + +import ( + "math" + "testing" +) + +// throwPanic will throw unexpected panics if true, or will just report them as errors if false +const throwPanic = true + +func dTolEqual(a, b float64) bool { + m := math.Max(math.Abs(a), math.Abs(b)) + if m > 1 { + a /= m + b /= m + } + if math.Abs(a-b) < 1e-14 { + return true + } + return false +} + +func dSliceTolEqual(a, b []float64) bool { + if len(a) != len(b) { + return false + } + for i := range a { + if !dTolEqual(a[i], b[i]) { + return false + } + } + return true +} + +func dSliceEqual(a, b []float64) bool { + if len(a) != len(b) { + return false + } + for i := range a { + if !(a[i] == b[i]) { + return false + } + } + return true +} + +func dCopyTwoTmp(x, xTmp, y, yTmp []float64) { + if len(x) != len(xTmp) { + panic("x size mismatch") + } + if len(y) != len(yTmp) { + panic("y size mismatch") + } + for i, val := range x { + xTmp[i] = val + } + for i, val := range y { + yTmp[i] = val + } +} + +// returns true if the function panics +func panics(f func()) (b bool) { + defer func() { + err := recover() + if err != nil { + b = true + } + }() + f() + return +} + +func testpanics(f func(), name string, t *testing.T) { + b := panics(f) + if !b { + t.Errorf("%v should panic and does not", name) + } +} diff --git a/testblas/level1double.go b/testblas/level1double.go index ca964fa4..4b62f7ed 100644 --- a/testblas/level1double.go +++ b/testblas/level1double.go @@ -9,76 +9,6 @@ import ( "testing" ) -func dTolEqual(a, b float64) bool { - m := math.Max(math.Abs(a), math.Abs(b)) - if m > 1 { - a /= m - b /= m - } - if math.Abs(a-b) < 1e-14 { - return true - } - return false -} - -func dSliceTolEqual(a, b []float64) bool { - if len(a) != len(b) { - return false - } - for i := range a { - if !dTolEqual(a[i], b[i]) { - return false - } - } - return true -} - -func dSliceEqual(a, b []float64) bool { - if len(a) != len(b) { - return false - } - for i := range a { - if !(a[i] == b[i]) { - return false - } - } - return true -} - -func dCopyTwoTmp(x, xTmp, y, yTmp []float64) { - if len(x) != len(xTmp) { - panic("x size mismatch") - } - if len(y) != len(yTmp) { - panic("y size mismatch") - } - for i, val := range x { - xTmp[i] = val - } - for i, val := range y { - yTmp[i] = val - } -} - -// returns true if the function panics -func panics(f func()) (b bool) { - defer func() { - err := recover() - if err != nil { - b = true - } - }() - f() - return -} - -func testpanics(f func(), name string, t *testing.T) { - b := panics(f) - if !b { - t.Errorf("%v should panic and does not", name) - } -} - type DoubleOneVectorCase struct { Name string X []float64 diff --git a/testblas/level2double.go b/testblas/level2double.go index cf9e17a9..b7cc2cb9 100644 --- a/testblas/level2double.go +++ b/testblas/level2double.go @@ -6,9 +6,6 @@ import ( "github.com/gonum/blas" ) -// throwPanic will throw unexpected panics if true, or will just report them as errors if false -const throwPanic = true - type DgemvCase struct { Name string m int