Add dlanst and test

native/dlanst.go

@@ -0,0 +1,68 @@
+// Copyright ©2016 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 native
+
+import (
+	"math"
+
+	"github.com/gonum/lapack"
+)
+
+// Dlanst computes the specified norm of a symmetric tridiagonal matrix A.
+// The diagonal elements of A are stored in d and the off-diagonal elements
+// are stored in e.
+func (impl Implementation) Dlanst(norm lapack.MatrixNorm, n int, d, e []float64) float64 {
+	if len(d) < n {
+		panic(badD)
+	}
+	if len(e) < n-1 {
+		panic(badE)
+	}
+	if n <= 0 {
+		return 0
+	}
+	switch norm {
+	default:
+		panic(badNorm)
+	case lapack.MaxAbs:
+		anorm := math.Abs(d[n-1])
+		for i := 0; i < n-1; i++ {
+			sum := math.Abs(d[i])
+			if anorm < sum || math.IsNaN(sum) {
+				anorm = sum
+			}
+			sum = math.Abs(e[i])
+			if anorm < sum || math.IsNaN(sum) {
+				anorm = sum
+			}
+		}
+		return anorm
+	case lapack.MaxColumnSum, lapack.MaxRowSum:
+		if n == 1 {
+			return math.Abs(d[0])
+		}
+		anorm := math.Abs(d[0]) + math.Abs(e[0])
+		sum := math.Abs(e[n-2]) + math.Abs(d[n-1])
+		if anorm < sum || math.IsNaN(sum) {
+			anorm = sum
+		}
+		for i := 1; i < n-1; i++ {
+			sum := math.Abs(d[i]) + math.Abs(e[i]) + math.Abs(e[i-1])
+			if anorm < sum || math.IsNaN(sum) {
+				anorm = sum
+			}
+		}
+		return anorm
+	case lapack.NormFrob:
+		var scale float64
+		sum := 1.0
+		if n > 1 {
+			scale, sum = impl.Dlassq(n-1, e, 1, scale, sum)
+			sum = 2 * sum
+		}
+		scale, sum = impl.Dlassq(n, d, 1, scale, sum)
+		return scale * math.Sqrt(sum)
+	}
+}

native/lapack_test.go

@@ -84,6 +84,10 @@ func TestDlas2(t *testing.T) {
 	testlapack.Dlas2Test(t, impl)
 }
 
+func TestDlanst(t *testing.T) {
+	testlapack.DlanstTest(t, impl)
+}
+
 func TestDlansy(t *testing.T) {
 	testlapack.DlansyTest(t, impl)
 }

testlapack/dlanst.go

@@ -0,0 +1,52 @@
+// Copyright ©2016 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 (
+	"math"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/lapack"
+)
+
+type Dlanster interface {
+	Dlanst(norm lapack.MatrixNorm, n int, d, e []float64) float64
+	Dlanger
+}
+
+func DlanstTest(t *testing.T, impl Dlanster) {
+	for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxColumnSum, lapack.MaxRowSum, lapack.NormFrob} {
+		for _, n := range []int{1, 3, 10, 100} {
+			for cas := 0; cas < 100; cas++ {
+				d := make([]float64, n)
+				for i := range d {
+					d[i] = rand.NormFloat64()
+				}
+				e := make([]float64, n-1)
+				for i := range e {
+					e[i] = rand.NormFloat64()
+				}
+
+				m := n
+				lda := n
+				a := make([]float64, m*lda)
+				for i := 0; i < n; i++ {
+					a[i*lda+i] = d[i]
+				}
+				for i := 0; i < n-1; i++ {
+					a[i*lda+i+1] = e[i]
+					a[(i+1)*lda+i] = e[i]
+				}
+				work := make([]float64, n)
+				syNorm := impl.Dlanst(norm, n, d, e)
+				geNorm := impl.Dlange(norm, m, n, a, lda, work)
+				if math.Abs(syNorm-geNorm) > 1e-12 {
+					t.Errorf("Norm mismatch: norm = %v, cas = %v, n = %v. Want %v, got %v.", string(norm), cas, n, geNorm, syNorm)
+				}
+			}
+		}
+	}
+}