lapack/{gonum,testlapack}: add Dlangb and its test

lapack/gonum/dlangb.go

@@ -0,0 +1,98 @@
+// 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 n×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, 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:
+		panic(badNorm)
+	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 n == 0 {
+		return 0
+	}
+
+	switch {
+	case len(ab) < (n-1)*ldab+ncol:
+		panic(shortAB)
+	case len(work) < n && norm == lapack.MaxColumnSum:
+		panic(shortWork)
+	}
+
+	var value float64
+	switch norm {
+	case lapack.MaxAbs:
+		for i := 0; i < n; 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 < n; 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:
+		work = work[:n]
+		for j := range work {
+			work[j] = 0
+		}
+		for i := 0; i < n; 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
+			}
+		}
+	case lapack.Frobenius:
+		scale := 0.0
+		ssq := 1.0
+		for i := 0; i < n; 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
+}

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"

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)

lapack/lapack64/lapack64.go

@@ -466,6 +466,20 @@ 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 square, n×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 {
+	if a.Rows != a.Cols {
+		panic("lapack64: matrix is not square")
+	}
+	return gonum.Implementation{}.Dlangb(norm, a.Rows, a.KL, a.KU, a.Data, a.Stride, work)
+}
+
 // 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

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, n, kl, ku int, ab []float64, ldab int, work []float64) 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 _, 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, n, kl, ku, ldab)
+							}
+						}
+					}
+				}
+			}
+		})
+	}
+}
+
+func dlangbTest(t *testing.T, impl Dlangber, rnd *rand.Rand, norm lapack.MatrixNorm, n, kl, ku, ldab int) {
+	const tol = 1e-14
+
+	name := fmt.Sprintf("n=%v,kl=%v,ku=%v,ldab=%v", n, kl, ku, ldab)
+
+	// Generate a random band matrix.
+	ab := randomSlice(n*ldab, rnd)
+	// Sometimes put a NaN into the matrix.
+	if n > 0 && rnd.Float64() < 0.5 {
+		i := rnd.Intn(n)
+		ab[i*ldab+kl] = math.NaN()
+	}
+	abCopy := make([]float64, len(ab))
+	copy(abCopy, ab)
+
+	// Deal with zero-sized matrices early.
+	if n == 0 {
+		got := impl.Dlangb(norm, n, kl, ku, nil, ldab, nil)
+		if got != 0 {
+			t.Errorf("%v: unexpected result for zero-sized matrix with nil input", name)
+		}
+		got = impl.Dlangb(norm, n, kl, ku, ab, ldab, nil)
+		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
+	}
+
+	// Generate a dense representation of the matrix and compute the wanted result.
+	a := zeros(n, n, n)
+	for i := 0; i < n; 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]
+		}
+	}
+
+	var work []float64
+	if norm == lapack.MaxColumnSum {
+		work = make([]float64, n)
+	}
+	got := impl.Dlangb(norm, n, kl, ku, ab, ldab, work)
+
+	if !floats.Same(ab, abCopy) {
+		t.Errorf("%v: unexpected modification in ab", name)
+	}
+
+	want := dlange(norm, n, n, 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 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)
+	}
+}