Repository: gonum/gonum
Pull request #1620 (2021-05-06 17:32:12)

lapack/{gonum,testlapack}: add Dlangb and its test
Preparation for adding `mat.Normer` and using it in `mat.Norm` similarly to `mat.Trace`.

Please take a look.

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->


codecov-commenter: # [Codecov](https://codecov.io/gh/gonum/gonum/pull/1620?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) Report
> Merging [#1620](https://codecov.io/gh/gonum/gonum/pull/1620?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (0b80ad8) into [master](https://codecov.io/gh/gonum/gonum/commit/3f0deff25eb9ebdd8a74ce9bd2754ec9c4b67a6e?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (3f0deff) will **decrease** coverage by `0.04%`.
> The diff coverage is `73.58%`.

[![Impacted file tree graph](https://codecov.io/gh/gonum/gonum/pull/1620/graphs/tree.svg?width=650&height=150&src=pr&token=B7qHnqSPNL&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)](https://codecov.io/gh/gonum/gonum/pull/1620?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)

```diff
@@            Coverage Diff             @@
##           master    #1620      +/-   ##
==========================================
- Coverage   73.20%   73.16%   -0.05%     
==========================================
  Files         502      498       -4     
  Lines       59510    59306     -204     
==========================================
- Hits        43563    43389     -174     
+ Misses      13491    13484       -7     
+ Partials     2456     2433      -23     
```


| [Impacted Files](https://codecov.io/gh/gonum/gonum/pull/1620?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | Coverage Δ | |
|---|---|---|
| [lapack/gonum/dlangb.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bGFwYWNrL2dvbnVtL2RsYW5nYi5nbw==) | `73.58% <73.58%> (ø)` | |
| [lapack/gonum/dlaqr1.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bGFwYWNrL2dvbnVtL2RsYXFyMS5nbw==) | `40.62% <0.00%> (-15.63%)` | :arrow_down: |
| [optimize/listsearch.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-b3B0aW1pemUvbGlzdHNlYXJjaC5nbw==) | `72.22% <0.00%> (-1.86%)` | :arrow_down: |
| [mat/index\_bound\_checks.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L2luZGV4X2JvdW5kX2NoZWNrcy5nbw==) | | |
| [graph/iterator/nodes\_map\_safe.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvaXRlcmF0b3Ivbm9kZXNfbWFwX3NhZmUuZ28=) | | |
| [...sp/fourier/internal/fftpack/array\_bounds\_checks.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-ZHNwL2ZvdXJpZXIvaW50ZXJuYWwvZmZ0cGFjay9hcnJheV9ib3VuZHNfY2hlY2tzLmdv) | | |
| [graph/internal/set/same\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvaW50ZXJuYWwvc2V0L3NhbWVfYXBwZW5naW5lLmdv) | | |
| [mat/offset\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L29mZnNldF9hcHBlbmdpbmUuZ28=) | | |
| [graph/graphs/gen/gen.go](https://codecov.io/gh/gonum/gonum/pull/1620/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvZ3JhcGhzL2dlbi9nZW4uZ28=) | `100.00% <0.00%> (+16.94%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gonum/gonum/pull/1620?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gonum/gonum/pull/1620?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Last update [3f0deff...0b80ad8](https://codecov.io/gh/gonum/gonum/pull/1620?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).

vladimir-ch: I don't understand how the codecov failure is happening.
kortschak: > I don't understand how the codecov failure is happening.

The parameter check panics are never exercised. This can be ignored.
vladimir-ch: > The parameter check panics are never exercised.

Ok, I misunderstood the codecov report. I thought it failed because of the negative delta in optimize.
kortschak: The change in optimize is not big enough to make a failure. It's irritating, because it's not consistent, but I haven't been able to stabilise it by repetition.

Patch:
diff --git a/lapack/gonum/dlangb.go b/lapack/gonum/dlangb.go
new file mode 100644
index 00000000..90e6efd0
--- /dev/null
+++ b/lapack/gonum/dlangb.go
@@ -0,0 +1,88 @@
+// Copyright ©2021 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 gonum
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/internal/asm/f64"
+	"gonum.org/v1/gonum/lapack"
+)
+
+// 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 {
+	ncol := kl + 1 + ku
+	switch {
+	case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
+		panic(badNorm)
+	case m < 0:
+		panic(mLT0)
+	case n < 0:
+		panic(nLT0)
+	case kl < 0:
+		panic(klLT0)
+	case ku < 0:
+		panic(kuLT0)
+	case ldab < ncol:
+		panic(badLdA)
+	}
+
+	// Quick return if possible.
+	if m == 0 || n == 0 {
+		return 0
+	}
+
+	switch {
+	case len(ab) < min(m, n+kl)*ldab:
+		panic(shortAB)
+	}
+
+	var value float64
+	switch norm {
+	case lapack.MaxAbs:
+		for i := 0; i < min(m, n+kl); i++ {
+			l := max(0, kl-i)
+			u := min(n+kl-i, ncol)
+			for _, aij := range ab[i*ldab+l : i*ldab+u] {
+				aij = math.Abs(aij)
+				if aij > value || math.IsNaN(aij) {
+					value = aij
+				}
+			}
+		}
+	case lapack.MaxRowSum:
+		for i := 0; i < min(m, n+kl); i++ {
+			l := max(0, kl-i)
+			u := min(n+kl-i, ncol)
+			sum := f64.L1Norm(ab[i*ldab+l : i*ldab+u])
+			if sum > value || math.IsNaN(sum) {
+				value = sum
+			}
+		}
+	case lapack.MaxColumnSum:
+		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:
+		scale := 0.0
+		ssq := 1.0
+		for i := 0; i < min(m, n+kl); i++ {
+			l := max(0, kl-i)
+			u := min(n+kl-i, ncol)
+			ilen := u - l
+			rowscale, rowssq := impl.Dlassq(ilen, ab[i*ldab+l:], 1, 0, 1)
+			scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+		}
+		value = scale * math.Sqrt(ssq)
+	}
+	return value
+}
diff --git a/lapack/gonum/errors.go b/lapack/gonum/errors.go
index 9d31e0d1..28f3b150 100644
--- a/lapack/gonum/errors.go
+++ b/lapack/gonum/errors.go
@@ -65,6 +65,8 @@ const (
 	kLT0        = "lapack: k < 0"
 	kLT1        = "lapack: k < 1"
 	kdLT0       = "lapack: kd < 0"
+	klLT0       = "lapack: kl < 0"
+	kuLT0       = "lapack: ku < 0"
 	mGTN        = "lapack: m > n"
 	mLT0        = "lapack: m < 0"
 	mmLT0       = "lapack: mm < 0"
diff --git a/lapack/lapack64/lapack64.go b/lapack/lapack64/lapack64.go
index c441bc7c..10724d49 100644
--- a/lapack/lapack64/lapack64.go
+++ b/lapack/lapack64/lapack64.go
@@ -466,6 +466,15 @@ func Lange(norm lapack.MatrixNorm, a blas64.General, work []float64) float64 {
 	return lapack64.Dlange(norm, a.Rows, a.Cols, a.Data, max(1, a.Stride), work)
 }
 
+// Langb returns the given norm of a general m×n band matrix with kl sub-diagonals and
+// ku super-diagonals.
+//
+// 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) 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.
 //
 // Dlangt is not part of the lapack.Float64 interface and so calls to Langt are


Test patch:
diff --git a/internal/testrand/extreme.go b/internal/testrand/extreme.go
index 97979c9e..3f4d3999 100644
--- a/internal/testrand/extreme.go
+++ b/internal/testrand/extreme.go
@@ -86,7 +86,7 @@ func (e *extreme) Float64() float64 {
 		return math.NaN()
 	}
 
-	return e.Float64()
+	return e.rnd.Float64()
 }
 
 // Int returns a non-negative pseudo-random int.
@@ -94,7 +94,7 @@ func (e *extreme) Int() int {
 	if e.p() {
 		return extremeInt[e.rnd.Intn(len(extremeInt))]
 	}
-	return e.Int()
+	return e.rnd.Int()
 }
 
 // Int31 returns a non-negative pseudo-random int32.
@@ -102,7 +102,7 @@ func (e *extreme) Int31() int32 {
 	if e.p() {
 		return extremeInt31[e.rnd.Intn(len(extremeInt31))]
 	}
-	return e.Int31()
+	return e.rnd.Int31()
 }
 
 // Int31n returns a non-negative pseudo-random int32 from range [0, n).
@@ -119,7 +119,7 @@ func (e *extreme) Int31n(n int32) int32 {
 			return n - 1
 		}
 	}
-	return e.Int31n(n)
+	return e.rnd.Int31n(n)
 }
 
 // Int63 returns a non-negative pseudo-random int64.
@@ -127,7 +127,7 @@ func (e *extreme) Int63() int64 {
 	if e.p() {
 		return extremeInt63[e.rnd.Intn(len(extremeInt63))]
 	}
-	return e.Int63()
+	return e.rnd.Int63()
 }
 
 // Int63n returns a non-negative pseudo-random int from range [0, n).
@@ -144,7 +144,7 @@ func (e *extreme) Int63n(n int64) int64 {
 			return n - 1
 		}
 	}
-	return e.Int63n(n)
+	return e.rnd.Int63n(n)
 }
 
 // Int returns a non-negative pseudo-random int from range [0, n).
@@ -161,7 +161,7 @@ func (e *extreme) Intn(n int) int {
 			return n - 1
 		}
 	}
-	return e.Intn(n)
+	return e.rnd.Intn(n)
 }
 
 // NormFloat64 returns a normally distributed pseudo-random float64 in range [-math.MaxFloat64, math.MaxFloat64].
@@ -173,7 +173,7 @@ func (e *extreme) NormFloat64() float64 {
 		return math.NaN()
 	}
 
-	return e.NormFloat64()
+	return e.rnd.NormFloat64()
 }
 
 // Uint32 returns a pseudo-random uint32.
@@ -181,7 +181,7 @@ func (e *extreme) Uint32() uint32 {
 	if e.p() {
 		return extremeUint32[e.rnd.Intn(len(extremeUint32))]
 	}
-	return e.Uint32()
+	return e.rnd.Uint32()
 }
 
 // Uint32 returns a pseudo-random uint64.
@@ -189,7 +189,7 @@ func (e *extreme) Uint64() uint64 {
 	if e.p() {
 		return extremeUint64[e.rnd.Intn(len(extremeUint64))]
 	}
-	return e.Uint64()
+	return e.rnd.Uint64()
 }
 
 // Uint64n returns a pseudo-random uint64 from range [0, n).
@@ -206,5 +206,5 @@ func (e *extreme) Uint64n(n uint64) uint64 {
 			return n - 1
 		}
 	}
-	return e.Uint64n(n)
+	return e.rnd.Uint64n(n)
 }
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index 4766d54a..4d76525b 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -208,6 +208,11 @@ func TestDlaln2(t *testing.T) {
 	testlapack.Dlaln2Test(t, impl)
 }
 
+func TestDlangb(t *testing.T) {
+	t.Parallel()
+	testlapack.DlangbTest(t, impl)
+}
+
 func TestDlange(t *testing.T) {
 	t.Parallel()
 	testlapack.DlangeTest(t, impl)
diff --git a/lapack/testlapack/dlangb.go b/lapack/testlapack/dlangb.go
new file mode 100644
index 00000000..574b785d
--- /dev/null
+++ b/lapack/testlapack/dlangb.go
@@ -0,0 +1,109 @@
+// Copyright ©2021 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 testlapack
+
+import (
+	"fmt"
+	"math"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/lapack"
+)
+
+type Dlangber interface {
+	Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab []float64, ldab int) float64
+}
+
+func DlangbTest(t *testing.T, impl Dlangber) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxRowSum, lapack.MaxColumnSum, lapack.Frobenius} {
+		t.Run(normToString(norm), func(t *testing.T) {
+			for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
+				for _, m := range []int{0, 1, 2, 3, 4, 5, 10} {
+					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} {
+								dlangbTest(t, impl, rnd, norm, m, n, kl, ku, ldab)
+							}
+						}
+					}
+				}
+			}
+		})
+	}
+}
+
+func dlangbTest(t *testing.T, impl Dlangber, rnd *rand.Rand, norm lapack.MatrixNorm, m, n, kl, ku, ldab int) {
+	const tol = 1e-14
+
+	name := fmt.Sprintf("m=%v,n=%v,kl=%v,ku=%v,ldab=%v", m, n, kl, ku, ldab)
+
+	// Generate a random band matrix.
+	ab := randomSlice(m*ldab, rnd)
+	// Sometimes put a NaN into the matrix.
+	if m > 0 && n > 0 && rnd.Float64() < 0.5 {
+		i := rnd.Intn(m)
+		ab[i*ldab+kl] = math.NaN()
+	}
+	abCopy := make([]float64, len(ab))
+	copy(abCopy, ab)
+
+	// Deal with zero-sized matrices early.
+	if m == 0 || n == 0 {
+		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)
+		if !floats.Same(ab, abCopy) {
+			t.Errorf("%v: unexpected modification in dl", name)
+		}
+		if got != 0 {
+			t.Errorf("%v: unexpected result for zero-sized matrix with non-nil input", name)
+		}
+		return
+	}
+
+	got := impl.Dlangb(norm, m, n, kl, ku, ab, ldab)
+
+	if !floats.Same(ab, abCopy) {
+		t.Errorf("%v: unexpected modification in ab", name)
+	}
+
+	// Generate a dense representation of the matrix and compute the wanted result.
+	a := zeros(m, n, n)
+	for i := 0; i < m; i++ {
+		for j := max(0, i-kl); j < min(i+ku+1, n); j++ {
+			a.Data[i*a.Stride+j] = ab[i*ldab+j-i+kl]
+		}
+	}
+	want := dlange(norm, a.Rows, a.Cols, a.Data, a.Stride)
+
+	if math.IsNaN(want) {
+		if !math.IsNaN(got) {
+			t.Errorf("%v: unexpected result with NaN element; got %v, want %v", name, got, want)
+		}
+		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)
+		}
+		return
+	}
+	diff := math.Abs(got - want)
+	if diff > tol {
+		t.Errorf("%v: unexpected result; got %v, want %v, diff=%v", name, got, want, diff)
+	}
+}