lapack/{gonum,testlapack}: don't use work in Dlangb

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.
-//
+func (impl Implementation) Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab []float64, ldab int) float64 {
-// 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 {
 	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 := 0; j < min(m+ku, n); j++ {
-		for j := range work {
+			jb := min(kl+j, ncol-1)
-			work[j] = 0
+			ib := max(0, j-ku)
-		}
+			jlen := min(j+kl, m-1) - ib + 1
-		for i := 0; i < min(m, n+kl); i++ {
+			sum := f64.L1NormInc(ab[ib*ldab+jb:], jlen, max(1, ldab-1))
-			l := max(0, kl-i)
+			if sum > value || math.IsNaN(sum) {
-			u := min(n+kl-i, ncol)
+				value = sum
-			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
 			}
 		}
 	case lapack.Frobenius:

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 {
+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), work)
+	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.

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,14 +28,12 @@ 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)
 							}
 						}
 					}
 				}
 			}
-			}
 		})
 	}
 }
@@ -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
+	got := impl.Dlangb(norm, m, n, kl, ku, ab, ldab)
-	if norm == lapack.MaxColumnSum {
-		work = make([]float64, n)
-	}
-	got := impl.Dlangb(norm, m, n, kl, ku, ab, ldab, work)
 
 	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)