From f0d77794af458e32971eb2aaa7f648331347ae7e Mon Sep 17 00:00:00 2001 From: Vladimir Chalupecky Date: Fri, 7 May 2021 11:57:25 +0200 Subject: [PATCH] lapack/{gonum,testlapack}: don't use work in Dlangb --- lapack/gonum/dlangb.go | 28 ++++++++-------------------- lapack/lapack64/lapack64.go | 6 ++---- lapack/testlapack/dlangb.go | 21 ++++++++++----------- 3 files changed, 20 insertions(+), 35 deletions(-) diff --git a/lapack/gonum/dlangb.go b/lapack/gonum/dlangb.go index 5bb8dfda..90e6efd0 100644 --- a/lapack/gonum/dlangb.go +++ b/lapack/gonum/dlangb.go @@ -13,9 +13,7 @@ import ( // Dlangb returns the given norm of an m×n band matrix with kl sub-diagonals and // ku super-diagonals. -// -// When norm is lapack.MaxColumnSum, the length of work must be at least n. -func (impl Implementation) Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab []float64, ldab int, work []float64) float64 { +func (impl Implementation) Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab []float64, ldab int) float64 { ncol := kl + 1 + ku switch { case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius: @@ -40,8 +38,6 @@ func (impl Implementation) Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab [ switch { case len(ab) < min(m, n+kl)*ldab: panic(shortAB) - case len(work) < n && norm == lapack.MaxColumnSum: - panic(shortWork) } var value float64 @@ -67,21 +63,13 @@ func (impl Implementation) Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab [ } } case lapack.MaxColumnSum: - work = work[:n] - for j := range work { - work[j] = 0 - } - for i := 0; i < min(m, n+kl); i++ { - l := max(0, kl-i) - u := min(n+kl-i, ncol) - for jb, aij := range ab[i*ldab+l : i*ldab+u] { - j := l + jb - kl + i - work[j] += math.Abs(aij) - } - } - for _, sumj := range work { - if sumj > value || math.IsNaN(sumj) { - value = sumj + for j := 0; j < min(m+ku, n); j++ { + jb := min(kl+j, ncol-1) + ib := max(0, j-ku) + jlen := min(j+kl, m-1) - ib + 1 + sum := f64.L1NormInc(ab[ib*ldab+jb:], jlen, max(1, ldab-1)) + if sum > value || math.IsNaN(sum) { + value = sum } } case lapack.Frobenius: diff --git a/lapack/lapack64/lapack64.go b/lapack/lapack64/lapack64.go index 990ce107..10724d49 100644 --- a/lapack/lapack64/lapack64.go +++ b/lapack/lapack64/lapack64.go @@ -469,12 +469,10 @@ func Lange(norm lapack.MatrixNorm, a blas64.General, work []float64) float64 { // Langb returns the given norm of a general m×n band matrix with kl sub-diagonals and // ku super-diagonals. // -// When norm is lapack.MaxColumnSum, the length of work must be at least n. -// // Dlangb is not part of the lapack.Float64 interface and so calls to Langb are always // executed by the Gonum implementation. -func Langb(norm lapack.MatrixNorm, a blas64.Band, work []float64) float64 { - return gonum.Implementation{}.Dlangb(norm, a.Rows, a.Cols, a.KL, a.KU, a.Data, max(1, a.Stride), work) +func Langb(norm lapack.MatrixNorm, a blas64.Band) float64 { + return gonum.Implementation{}.Dlangb(norm, a.Rows, a.Cols, a.KL, a.KU, a.Data, max(1, a.Stride)) } // Langt computes the specified norm of an n×n tridiagonal matrix. diff --git a/lapack/testlapack/dlangb.go b/lapack/testlapack/dlangb.go index 77c9cf98..574b785d 100644 --- a/lapack/testlapack/dlangb.go +++ b/lapack/testlapack/dlangb.go @@ -16,7 +16,7 @@ import ( ) type Dlangber interface { - Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab []float64, ldab int, work []float64) float64 + Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab []float64, ldab int) float64 } func DlangbTest(t *testing.T, impl Dlangber) { @@ -28,9 +28,7 @@ func DlangbTest(t *testing.T, impl Dlangber) { for _, kl := range []int{0, 1, 2, 3, 4, 5, 10} { for _, ku := range []int{0, 1, 2, 3, 4, 5, 10} { for _, ldab := range []int{kl + ku + 1, kl + ku + 1 + 7} { - for iter := 0; iter < 10; iter++ { - dlangbTest(t, impl, rnd, norm, m, n, kl, ku, ldab) - } + dlangbTest(t, impl, rnd, norm, m, n, kl, ku, ldab) } } } @@ -57,11 +55,11 @@ func dlangbTest(t *testing.T, impl Dlangber, rnd *rand.Rand, norm lapack.MatrixN // Deal with zero-sized matrices early. if m == 0 || n == 0 { - got := impl.Dlangb(norm, m, n, kl, ku, nil, ldab, nil) + got := impl.Dlangb(norm, m, n, kl, ku, nil, ldab) if got != 0 { t.Errorf("%v: unexpected result for zero-sized matrix with nil input", name) } - got = impl.Dlangb(norm, m, n, kl, ku, ab, ldab, nil) + got = impl.Dlangb(norm, m, n, kl, ku, ab, ldab) if !floats.Same(ab, abCopy) { t.Errorf("%v: unexpected modification in dl", name) } @@ -71,11 +69,7 @@ func dlangbTest(t *testing.T, impl Dlangber, rnd *rand.Rand, norm lapack.MatrixN return } - var work []float64 - if norm == lapack.MaxColumnSum { - work = make([]float64, n) - } - got := impl.Dlangb(norm, m, n, kl, ku, ab, ldab, work) + got := impl.Dlangb(norm, m, n, kl, ku, ab, ldab) if !floats.Same(ab, abCopy) { t.Errorf("%v: unexpected modification in ab", name) @@ -97,6 +91,11 @@ func dlangbTest(t *testing.T, impl Dlangber, rnd *rand.Rand, norm lapack.MatrixN return } + if math.IsNaN(got) { + t.Errorf("%v: unexpected NaN; want %v", name, want) + return + } + if norm == lapack.MaxAbs { if got != want { t.Errorf("%v: unexpected result; got %v, want %v", name, got, want)