diff --git a/native/dlanst.go b/native/dlanst.go new file mode 100644 index 0000000..d806de9 --- /dev/null +++ b/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) + } +} diff --git a/native/lapack_test.go b/native/lapack_test.go index 60387b3..424a837 100644 --- a/native/lapack_test.go +++ b/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) } diff --git a/testlapack/dlanst.go b/testlapack/dlanst.go new file mode 100644 index 0000000..264ec68 --- /dev/null +++ b/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) + } + } + } + } +}