lapack: add Dlangt

lapack/gonum/dlangt.go

@@ -0,0 +1,115 @@
+// Copyright ©2020 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/lapack"
+)
+
+// Dlangt returns the value of the given norm of an n×n tridiagonal matrix
+// represented by the three diagonals.
+//
+// d must have length at least n and dl and du must have length at least n-1.
+func (impl Implementation) Dlangt(norm lapack.MatrixNorm, n int, dl, d, du []float64) float64 {
+	switch {
+	case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
+		panic(badNorm)
+	case n < 0:
+		panic(nLT0)
+	}
+
+	if n == 0 {
+		return 0
+	}
+
+	switch {
+	case len(dl) < n-1:
+		panic(shortDL)
+	case len(d) < n:
+		panic(shortD)
+	case len(du) < n-1:
+		panic(shortDU)
+	}
+
+	dl = dl[:n-1]
+	d = d[:n]
+	du = du[:n-1]
+
+	var anorm float64
+	switch norm {
+	case lapack.MaxAbs:
+		for _, diag := range [][]float64{dl, d, du} {
+			for _, di := range diag {
+				if math.IsNaN(di) {
+					return di
+				}
+				di = math.Abs(di)
+				if di > anorm {
+					anorm = di
+				}
+			}
+		}
+	case lapack.MaxColumnSum:
+		if n == 1 {
+			return math.Abs(d[0])
+		}
+		anorm = math.Abs(d[0]) + math.Abs(dl[0])
+		if math.IsNaN(anorm) {
+			return anorm
+		}
+		tmp := math.Abs(du[n-2]) + math.Abs(d[n-1])
+		if math.IsNaN(tmp) {
+			return tmp
+		}
+		if tmp > anorm {
+			anorm = tmp
+		}
+		for i := 1; i < n-1; i++ {
+			tmp = math.Abs(du[i-1]) + math.Abs(d[i]) + math.Abs(dl[i])
+			if math.IsNaN(tmp) {
+				return tmp
+			}
+			if tmp > anorm {
+				anorm = tmp
+			}
+		}
+	case lapack.MaxRowSum:
+		if n == 1 {
+			return math.Abs(d[0])
+		}
+		anorm = math.Abs(d[0]) + math.Abs(du[0])
+		if math.IsNaN(anorm) {
+			return anorm
+		}
+		tmp := math.Abs(dl[n-2]) + math.Abs(d[n-1])
+		if math.IsNaN(tmp) {
+			return tmp
+		}
+		if tmp > anorm {
+			anorm = tmp
+		}
+		for i := 1; i < n-1; i++ {
+			tmp = math.Abs(dl[i-1]) + math.Abs(d[i]) + math.Abs(du[i])
+			if math.IsNaN(tmp) {
+				return tmp
+			}
+			if tmp > anorm {
+				anorm = tmp
+			}
+		}
+	case lapack.Frobenius:
+		scale := 0.0
+		ssq := 1.0
+		scale, ssq = impl.Dlassq(n, d, 1, scale, ssq)
+		if n > 1 {
+			scale, ssq = impl.Dlassq(n-1, dl, 1, scale, ssq)
+			scale, ssq = impl.Dlassq(n-1, du, 1, scale, ssq)
+		}
+		anorm = scale * math.Sqrt(ssq)
+	}
+	return anorm
+}

lapack/gonum/errors.go

@@ -116,6 +116,8 @@ const (
 	shortC     = "lapack: insufficient length of c"
 	shortCNorm = "lapack: insufficient length of cnorm"
 	shortD     = "lapack: insufficient length of d"
+	shortDL    = "lapack: insufficient length of dl"
+	shortDU    = "lapack: insufficient length of du"
 	shortE     = "lapack: insufficient length of e"
 	shortF     = "lapack: insufficient length of f"
 	shortH     = "lapack: insufficient length of h"

lapack/gonum/lapack_test.go

@@ -203,6 +203,11 @@ func TestDlange(t *testing.T) {
 	testlapack.DlangeTest(t, impl)
 }
 
+func TestDlangt(t *testing.T) {
+	t.Parallel()
+	testlapack.DlangtTest(t, impl)
+}
+
 func TestDlapy2(t *testing.T) {
 	t.Parallel()
 	testlapack.Dlapy2Test(t, impl)

lapack/lapack64/lapack64.go

@@ -432,6 +432,14 @@ 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)
 }
 
+// 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
+// always executed by the Gonum implementation.
+func Langt(norm lapack.MatrixNorm, a Tridiagonal) float64 {
+	return gonum.Implementation{}.Dlangt(norm, a.N, a.DL, a.D, a.DU)
+}
+
 // Lansb computes the specified norm of an n×n symmetric band matrix. If
 // norm == lapack.MaxColumnSum or norm == lapack.MaxRowSum, work must have length
 // at least n and this function will panic otherwise.

lapack/testlapack/dlangt.go

@@ -0,0 +1,122 @@
+// Copyright ©2020 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 Dlangter interface {
+	Dlangt(norm lapack.MatrixNorm, n int, dl, d, du []float64) float64
+}
+
+func DlangtTest(t *testing.T, impl Dlangter) {
+	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 iter := 0; iter < 10; iter++ {
+					dlangtTest(t, impl, rnd, norm, n)
+				}
+			}
+		})
+	}
+}
+
+func dlangtTest(t *testing.T, impl Dlangter, rnd *rand.Rand, norm lapack.MatrixNorm, n int) {
+	const (
+		tol   = 1e-14
+		extra = 10
+	)
+
+	name := fmt.Sprintf("n=%v", n)
+
+	// Generate three random diagonals.
+	dl := randomSlice(n+extra, rnd)
+	dlCopy := make([]float64, len(dl))
+	copy(dlCopy, dl)
+
+	d := randomSlice(n+1+extra, rnd)
+	// Sometimes put a NaN into the matrix.
+	if n > 0 && rnd.Float64() < 0.5 {
+		d[rnd.Intn(n)] = math.NaN()
+	}
+	dCopy := make([]float64, len(d))
+	copy(dCopy, d)
+
+	du := randomSlice(n+extra, rnd)
+	duCopy := make([]float64, len(du))
+	copy(duCopy, du)
+
+	// Deal with zero-sized matrices early.
+	if n == 0 {
+		got := impl.Dlangt(norm, n, nil, nil, nil)
+		if got != 0 {
+			t.Errorf("%v: unexpected result for zero-sized matrix with nil input", name)
+		}
+		got = impl.Dlangt(norm, n, dl, d, du)
+		if !floats.Same(dl, dlCopy) {
+			t.Errorf("%v: unexpected modification in dl", name)
+		}
+		if !floats.Same(d, dCopy) {
+			t.Errorf("%v: unexpected modification in d", name)
+		}
+		if !floats.Same(du, duCopy) {
+			t.Errorf("%v: unexpected modification in du", 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-1; i++ {
+		a.Data[i*a.Stride+i] = d[i]
+		a.Data[i*a.Stride+i+1] = du[i]
+		a.Data[(i+1)*a.Stride+i] = dl[i]
+	}
+	a.Data[(n-1)*a.Stride+n-1] = d[n-1]
+
+	got := impl.Dlangt(norm, n, dl, d, du)
+
+	if !floats.Same(dl, dlCopy) {
+		t.Errorf("%v: unexpected modification in dl", name)
+	}
+	if !floats.Same(d, dCopy) {
+		t.Errorf("%v: unexpected modification in d", name)
+	}
+	if !floats.Same(du, duCopy) {
+		t.Errorf("%v: unexpected modification in du", 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)
+	}
+}