Repository: gonum/gonum
Pull request #1887 (2023-08-30 22:14:08)

lapack: add `Dgghrd`
Part of #1651. Test not added since it is not readily apparent to me how to go about testing these functions. My background is not in linear algebra. 

- Added `Dgghrd` with following signature:
   ```go
   func (Implementation) Dgghrd(compq, compz lapack.OrthoComp, n, ilo, ihi int, a []float64, lda int, b []float64, ldb int, q []float64, ldq int, z []float64, ldz int)
   ```
- Added `lapack.OrthoComp` type and `badOrthoComp` panic string.


codecov-commenter: ## [Codecov](https://app.codecov.io/gh/gonum/gonum/pull/1887?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) Report
Patch coverage has no change and project coverage change: **`-0.23%`** :warning:
> Comparison is base [(`93223b4`)](https://app.codecov.io/gh/gonum/gonum/commit/93223b4004117a4e189a5629056123931ce85314?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) 74.00% compared to head [(`c2449e3`)](https://app.codecov.io/gh/gonum/gonum/pull/1887?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) 73.78%.

:exclamation: Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please [install the GitHub App Integration](https://github.com/apps/codecov) for your organization. [Read more](https://about.codecov.io/blog/codecov-is-updating-its-github-integration/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #1887      +/-   ##
==========================================
- Coverage   74.00%   73.78%   -0.23%     
==========================================
  Files         526      509      -17     
  Lines       63210    62620     -590     
==========================================
- Hits        46778    46202     -576     
- Misses      13854    13873      +19     
+ Partials     2578     2545      -33     
```


| [Files Changed](https://app.codecov.io/gh/gonum/gonum/pull/1887?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | Coverage Δ | |
|---|---|---|
| [blas/blas64/blas64.go](https://app.codecov.io/gh/gonum/gonum/pull/1887?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Ymxhcy9ibGFzNjQvYmxhczY0Lmdv) | `72.50% <ø> (ø)` | |
| [lapack/gonum/dgghrd.go](https://app.codecov.io/gh/gonum/gonum/pull/1887?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bGFwYWNrL2dvbnVtL2RnZ2hyZC5nbw==) | `0.00% <0.00%> (ø)` | |

... and [21 files with indirect coverage changes](https://app.codecov.io/gh/gonum/gonum/pull/1887/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gonum/gonum/pull/1887?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).

soypat: > This needs a test.

Added test. did not test for ilo,ihi combo that are best suited for use with Dggbal.
soypat: Are there any blockers left in this PR?

Patch:
diff --git a/AUTHORS b/AUTHORS
index 91467eaa..8f3e4cc2 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -65,6 +65,7 @@ Jonathan J Lawlor <jonathan.lawlor@gmail.com>
 Jonathan Reiter <jonreiter@gmail.com>
 Jonathan Schroeder <jd.schroeder@gmail.com>
 Joost van Amersfoort <git@joo.st>
+Jordan Stoker <jordan_stoker@hotmail.com>
 Joseph Watson <jtwatson@linux-consulting.us>
 Josh Wilson <josh.craig.wilson@gmail.com>
 Julien Roland <juroland@gmail.com>
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index 93e01a2c..b7d499d6 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -73,6 +73,7 @@ Jonathan J Lawlor <jonathan.lawlor@gmail.com>
 Jonathan Reiter <jonreiter@gmail.com>
 Jonathan Schroeder <jd.schroeder@gmail.com>
 Joost van Amersfoort <git@joo.st>
+Jordan Stoker <jordan_stoker@hotmail.com>
 Joseph Watson <jtwatson@linux-consulting.us>
 Josh Wilson <josh.craig.wilson@gmail.com>
 Julien Roland <juroland@gmail.com>
diff --git a/blas/blas64/blas64.go b/blas/blas64/blas64.go
index c336dc89..64ac985c 100644
--- a/blas/blas64/blas64.go
+++ b/blas/blas64/blas64.go
@@ -20,7 +20,7 @@ func Use(b blas.Float64) {
 
 // Implementation returns the current BLAS float64 implementation.
 //
-// Implementation allows direct calls to the current the BLAS float64 implementation
+// Implementation allows direct calls to the current BLAS float64 implementation
 // giving finer control of parameters.
 func Implementation() blas.Float64 {
 	return blas64
diff --git a/lapack/gonum/dgghrd.go b/lapack/gonum/dgghrd.go
new file mode 100644
index 00000000..da09ea9e
--- /dev/null
+++ b/lapack/gonum/dgghrd.go
@@ -0,0 +1,125 @@
+// Copyright ©2023 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 (
+	"gonum.org/v1/gonum/blas"
+	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/lapack"
+)
+
+// Dgghrd reduces a pair of real matrices (A,B) to generalized upper
+// Hessenberg form using orthogonal transformations, where A is a
+// general matrix and B is upper triangular.  The form of the
+// generalized eigenvalue problem is
+//
+//	A*x = lambda*B*x,
+//
+// and B is typically made upper triangular by computing its QR
+// factorization and moving the orthogonal matrix Q to the left side
+// of the equation.
+// This subroutine simultaneously reduces A to a Hessenberg matrix H:
+//
+//	Qᵀ*A*Z = H
+//
+// and transforms B to another upper triangular matrix T:
+//
+//	Qᵀ*B*Z = T
+//
+// in order to reduce the problem to its standard form
+//
+//	H*y = lambda*T*y
+//
+// where y = Zᵀ*x.
+//
+// The orthogonal matrices Q and Z are determined as products of Givens
+// rotations.  They may either be formed explicitly, or they may be
+// postmultiplied into input matrices Q1 and Z1, so that
+//
+//	Q1 * A * Z1ᵀ = (Q1*Q) * H * (Z1*Z)ᵀ
+//	Q1 * B * Z1ᵀ = (Q1*Q) * T * (Z1*Z)ᵀ
+//
+// If Q1 is the orthogonal matrix from the QR factorization of B in the
+// original equation A*x = lambda*B*x, then Dgghrd reduces the original
+// problem to generalized Hessenberg form.
+//
+// Dgghrd is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dgghrd(compq, compz lapack.OrthoComp, n, ilo, ihi int, a []float64, lda int, b []float64, ldb int, q []float64, ldq int, z []float64, ldz int) {
+	switch {
+	case compq != lapack.OrthoNone && compq != lapack.OrthoEntry && compq != lapack.OrthoUnit:
+		panic(badOrthoComp)
+	case compz != lapack.OrthoNone && compz != lapack.OrthoEntry && compz != lapack.OrthoUnit:
+		panic(badOrthoComp)
+	case len(a) < (n-1)*lda+n:
+		panic(shortA)
+	case len(b) < (n-1)*ldb+n:
+		panic(shortB)
+	case n < 0:
+		panic(nLT0)
+	case ilo < 0:
+		panic(badIlo)
+	case ihi < ilo-1 || ihi >= n:
+		panic(badIhi)
+	case lda < max(1, n):
+		panic(badLdA)
+	case ldb < max(1, n):
+		panic(badLdB)
+	case (compq != lapack.OrthoNone && ldq < n) || ldq < 1:
+		panic(badLdQ)
+	case (compz != lapack.OrthoNone && ldz < n) || ldz < 1:
+		panic(badLdZ)
+	case compq != lapack.OrthoNone && len(q) < (n-1)*ldq+n:
+		panic(shortQ)
+	case compz != lapack.OrthoNone && len(z) < (n-1)*ldz+n:
+		panic(shortZ)
+	}
+
+	if compq == lapack.OrthoUnit {
+		impl.Dlaset(blas.All, n, n, 0, 1, q, ldq)
+	}
+	if compz == lapack.OrthoUnit {
+		impl.Dlaset(blas.All, n, n, 0, 1, z, ldz)
+	}
+	if n <= 1 {
+		return // Quick return if possible.
+	}
+
+	// Zero out lower triangle of B.
+	for i := 1; i < n; i++ {
+		for j := 0; j < i; j++ {
+			b[i*ldb+j] = 0
+		}
+	}
+	bi := blas64.Implementation()
+	// Reduce A and B.
+	for jcol := ilo; jcol <= ihi-2; jcol++ {
+		for jrow := ihi; jrow >= jcol+2; jrow-- {
+			// Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL).
+			var c, s float64
+			c, s, a[(jrow-1)*lda+jcol] = impl.Dlartg(a[(jrow-1)*lda+jcol], a[jrow*lda+jcol])
+			a[jrow*lda+jcol] = 0
+			bi.Drot(n-jcol-1, a[(jrow-1)*lda+jcol+1:], 1,
+				a[jrow*lda+jcol+1:], 1, c, s)
+
+			bi.Drot(n+2-jrow-1, b[(jrow-1)*ldb+jrow-1:], 1,
+				b[jrow*ldb+jrow-1:], 1, c, s)
+
+			if compq != lapack.OrthoNone {
+				bi.Drot(n, q[jrow-1:], ldq, q[jrow:], ldq, c, s)
+			}
+
+			// Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1).
+			c, s, b[jrow*ldb+jrow] = impl.Dlartg(b[jrow*ldb+jrow], b[jrow*ldb+jrow-1])
+			b[jrow*ldb+jrow-1] = 0
+
+			bi.Drot(ihi+1, a[jrow:], lda, a[jrow-1:], lda, c, s)
+			bi.Drot(jrow, b[jrow:], ldb, b[jrow-1:], ldb, c, s)
+
+			if compz != lapack.OrthoNone {
+				bi.Drot(n, z[jrow:], ldz, z[jrow-1:], ldz, c, s)
+			}
+		}
+	}
+}
diff --git a/lapack/gonum/dlanhs.go b/lapack/gonum/dlanhs.go
new file mode 100644
index 00000000..054b90f0
--- /dev/null
+++ b/lapack/gonum/dlanhs.go
@@ -0,0 +1,78 @@
+// Copyright ©2023 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/blas/blas64"
+	"gonum.org/v1/gonum/lapack"
+)
+
+// Dlanhs returns the value of the one norm, or the Frobenius norm, or the
+// infinity norm, or the element of largest absolute value of a Hessenberg
+// matrix A.
+//
+// If norm is lapack.MaxColumnSum, work must have length at least n.
+func (impl Implementation) Dlanhs(norm lapack.MatrixNorm, n int, a []float64, lda int, work []float64) float64 {
+	switch {
+	case norm != lapack.MaxRowSum && norm != lapack.MaxAbs && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
+		panic(badNorm)
+	case n < 0:
+		panic(nLT0)
+	case lda < max(1, n):
+		panic(badLdA)
+	}
+
+	if n == 0 {
+		return 0
+	}
+
+	switch {
+	case len(a) < (n-1)*lda+n:
+		panic(shortA)
+	case norm == lapack.MaxColumnSum && len(work) < n:
+		panic(shortWork)
+	}
+
+	bi := blas64.Implementation()
+	var value float64
+	switch norm {
+	case lapack.MaxAbs:
+		for i := 0; i < n; i++ {
+			minj := max(0, i-1)
+			for _, v := range a[i*lda+minj : i*lda+n] {
+				value = math.Max(value, math.Abs(v))
+			}
+		}
+	case lapack.MaxColumnSum:
+		for i := 0; i < n; i++ {
+			work[i] = 0
+		}
+		for i := 0; i < n; i++ {
+			for j := max(0, i-1); j < n; j++ {
+				work[j] += math.Abs(a[i*lda+j])
+			}
+		}
+		for _, v := range work[:n] {
+			value = math.Max(value, v)
+		}
+	case lapack.MaxRowSum:
+		for i := 0; i < n; i++ {
+			minj := max(0, i-1)
+			sum := bi.Dasum(n-minj, a[i*lda+minj:], 1)
+			value = math.Max(value, sum)
+		}
+	case lapack.Frobenius:
+		scale := 0.0
+		sum := 1.0
+		for i := 0; i < n; i++ {
+			minj := max(0, i-1)
+			scale, sum = impl.Dlassq(n-minj, a[i*lda+minj:], 1, scale, sum)
+		}
+		value = scale * math.Sqrt(sum)
+	}
+	return value
+}
diff --git a/lapack/gonum/errors.go b/lapack/gonum/errors.go
index 47652841..711cc2d5 100644
--- a/lapack/gonum/errors.go
+++ b/lapack/gonum/errors.go
@@ -21,6 +21,7 @@ const (
 	badMatrixType       = "lapack: bad MatrixType"
 	badMaximizeNormXJob = "lapack: bad MaximizeNormXJob"
 	badNorm             = "lapack: bad Norm"
+	badOrthoComp        = "lapack: bad OrthoComp"
 	badPivot            = "lapack: bad Pivot"
 	badRightEVJob       = "lapack: bad RightEVJob"
 	badSVDJob           = "lapack: bad SVDJob"
diff --git a/lapack/lapack.go b/lapack/lapack.go
index 5f60438f..72111d66 100644
--- a/lapack/lapack.go
+++ b/lapack/lapack.go
@@ -226,3 +226,12 @@ const (
 	LocalLookAhead       MaximizeNormXJob = 0 // Solve Z*x=h-f where h is a vector of ±1.
 	NormalizedNullVector MaximizeNormXJob = 2 // Compute an approximate null-vector e of Z, normalize e and solve Z*x=±e-f.
 )
+
+// OrthoComp specifies whether and how the orthogonal matrix is computed in Dgghrd.
+type OrthoComp byte
+
+const (
+	OrthoNone  OrthoComp = 'N' // Do not compute orthogonal matrix.
+	OrthoUnit  OrthoComp = 'I' // Argument is initialized to the unit matrix and the orthogonal matrix is returned.
+	OrthoEntry OrthoComp = 'V' // Argument Q contains orthogonal matrix Q1 on entry and the product Q1*Q is returned.
+)
diff --git a/spatial/kdtree/kdtree.go b/spatial/kdtree/kdtree.go
index 1bcfe1a3..0e1c522a 100644
--- a/spatial/kdtree/kdtree.go
+++ b/spatial/kdtree/kdtree.go
@@ -324,9 +324,13 @@ func NewNKeeper(n int) *NKeeper {
 func (k *NKeeper) Keep(c ComparableDist) {
 	if c.Dist <= k.Heap[0].Dist { // Favour later finds to displace sentinel.
 		if len(k.Heap) == cap(k.Heap) {
-			heap.Pop(k)
+			// Changing the value of the element at index 0 and then calling Fix is equivalent to,
+			// but less expensive than, calling Pop() followed by a Push of the new value.
+			k.Heap[0] = c
+			heap.Fix(k, 0)
+		} else {
+			heap.Push(k, c)
 		}
-		heap.Push(k, c)
 	}
 }
 


Test patch:
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index 5a3c418c..55d50616 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -148,6 +148,11 @@ func TestDgetrs(t *testing.T) {
 	testlapack.DgetrsTest(t, impl)
 }
 
+func TestDgghrd(t *testing.T) {
+	t.Parallel()
+	testlapack.DgghrdTest(t, impl)
+}
+
 func TestDggsvd3(t *testing.T) {
 	t.Parallel()
 	testlapack.Dggsvd3Test(t, impl)
@@ -238,6 +243,11 @@ func TestDlangt(t *testing.T) {
 	testlapack.DlangtTest(t, impl)
 }
 
+func TestDlanhs(t *testing.T) {
+	t.Parallel()
+	testlapack.DlanhsTest(t, impl)
+}
+
 func TestDlapy2(t *testing.T) {
 	t.Parallel()
 	testlapack.Dlapy2Test(t, impl)
diff --git a/lapack/testlapack/dgghrd.go b/lapack/testlapack/dgghrd.go
new file mode 100644
index 00000000..02b0bc86
--- /dev/null
+++ b/lapack/testlapack/dgghrd.go
@@ -0,0 +1,141 @@
+// Copyright ©2023 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"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/blas"
+	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/lapack"
+)
+
+type Dgghrder interface {
+	Dgghrd(compq, compz lapack.OrthoComp, n, ilo, ihi int, a []float64, lda int, b []float64, ldb int, q []float64, ldq int, z []float64, ldz int)
+}
+
+func DgghrdTest(t *testing.T, impl Dgghrder) {
+	const tol = 1e-13
+	const ldAdd = 5
+	rnd := rand.New(rand.NewSource(1))
+	comps := []lapack.OrthoComp{lapack.OrthoUnit, lapack.OrthoNone, lapack.OrthoEntry}
+	for _, compq := range comps {
+		for _, compz := range comps {
+			for _, n := range []int{2, 0, 1, 4, 15} {
+				ldMin := max(1, n)
+				for _, lda := range []int{ldMin, ldMin + ldAdd} {
+					for _, ldb := range []int{ldMin, ldMin + ldAdd} {
+						for _, ldq := range []int{ldMin, ldMin + ldAdd} {
+							for _, ldz := range []int{ldMin, ldMin + ldAdd} {
+								testDgghrd(t, impl, rnd, tol, compq, compz, n, 0, n-1, lda, ldb, ldq, ldz)
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
+func testDgghrd(t *testing.T, impl Dgghrder, rnd *rand.Rand, tol float64, compq, compz lapack.OrthoComp, n, ilo, ihi, lda, ldb, ldq, ldz int) {
+	a := randomGeneral(n, n, lda, rnd)
+	b := blockedUpperTriGeneral(n, n, 0, n, ldb, false, rnd)
+	var q, q1, z, z1 blas64.General
+	if compq == lapack.OrthoEntry {
+		q = randomOrthogonal(n, rnd)
+		q1 = cloneGeneral(q)
+	} else {
+		q = nanGeneral(n, n, ldq)
+	}
+	if compz == lapack.OrthoEntry {
+		z = randomOrthogonal(n, rnd)
+		z1 = cloneGeneral(z)
+	} else {
+		z = nanGeneral(n, n, ldz)
+	}
+
+	hGot := cloneGeneral(a)
+	tGot := cloneGeneral(b)
+	for i := 1; i < n; i++ {
+		for j := 0; j < i; j++ {
+			// Set all lower tri elems to NaN to catch bad implementations.
+			tGot.Data[i*tGot.Stride+j] = math.NaN()
+		}
+	}
+	impl.Dgghrd(compq, compz, n, ilo, ihi, hGot.Data, hGot.Stride, tGot.Data, tGot.Stride, q.Data, q.Stride, z.Data, z.Stride)
+	if n == 0 {
+		return
+	}
+	if !isUpperHessenberg(hGot) {
+		t.Error("H is not upper Hessenberg")
+	}
+	if !isNaNFree(tGot) || !isNaNFree(hGot) {
+		t.Error("T or H is/or not NaN free")
+	}
+	if !isUpperTriangular(tGot) {
+		t.Error("T is not upper triangular")
+	}
+	if compq == lapack.OrthoNone {
+		if !isAllNaN(q.Data) {
+			t.Errorf("Q is not NaN")
+		}
+		return
+	}
+	if compz == lapack.OrthoNone {
+		if !isAllNaN(z.Data) {
+			t.Errorf("Z is not NaN")
+		}
+		return
+	}
+	if compq != compz {
+		return // Do not handle mixed case
+	}
+	comp := compq
+	aux := zeros(n, n, n)
+
+	switch comp {
+	case lapack.OrthoUnit:
+		// Qᵀ*A*Z = H
+		hCalc := zeros(n, n, n)
+		blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, a, 0, aux)
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aux, z, 1, hCalc)
+		if !equalApproxGeneral(hGot, hCalc, tol) {
+			t.Errorf("Qᵀ*A*Z != H")
+		}
+
+		// Qᵀ*B*Z = T
+		tCalc := zeros(n, n, n)
+		blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, b, 0, aux)
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aux, z, 1, tCalc)
+		if !equalApproxGeneral(hGot, hCalc, tol) {
+			t.Errorf("Qᵀ*B*Z != T")
+		}
+	case lapack.OrthoEntry:
+		//	Q1 * A * Z1ᵀ = (Q1*Q) * H * (Z1*Z)ᵀ
+		lhs := zeros(n, n, n)
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q1, a, 0, aux)
+		blas64.Gemm(blas.NoTrans, blas.Trans, 1, aux, z1, 0, lhs) // lhs = Q1 * A * Z1ᵀ
+
+		rhs := zeros(n, n, n)
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, hGot, 0, aux)
+		blas64.Gemm(blas.NoTrans, blas.Trans, 1, aux, z, 0, rhs)
+		if !equalApproxGeneral(lhs, rhs, tol) {
+			t.Errorf("Q1 * A * Z1ᵀ != (Q1*Q) * H * (Z1*Z)ᵀ")
+		}
+
+		//	Q1 * B * Z1ᵀ = (Q1*Q) * T * (Z1*Z)ᵀ
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q1, b, 0, aux)
+		blas64.Gemm(blas.NoTrans, blas.Trans, 1, aux, z1, 0, lhs)
+
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, tGot, 0, aux)
+		blas64.Gemm(blas.NoTrans, blas.Trans, 1, aux, z, 0, rhs)
+		if !equalApproxGeneral(lhs, rhs, tol) {
+			t.Errorf("Q1 * B * Z1ᵀ != (Q1*Q) * T * (Z1*Z)ᵀ")
+		}
+	}
+}
diff --git a/lapack/testlapack/dlanhs.go b/lapack/testlapack/dlanhs.go
new file mode 100644
index 00000000..2c5bb2e3
--- /dev/null
+++ b/lapack/testlapack/dlanhs.go
@@ -0,0 +1,48 @@
+// Copyright ©2023 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"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/lapack"
+)
+
+type Dlanhser interface {
+	Dlanhs(norm lapack.MatrixNorm, n int, a []float64, lda int, work []float64) float64
+}
+
+func DlanhsTest(t *testing.T, impl Dlanhser) {
+	const tol = 1e-15
+	rnd := rand.New(rand.NewSource(1))
+	for _, n := range []int{0, 1, 2, 4, 9} {
+		for _, lda := range []int{max(1, n), n + 5} {
+			a := randomGeneral(n, n, lda, rnd)
+			for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxRowSum, lapack.MaxColumnSum, lapack.Frobenius} {
+				var work []float64
+				if norm == lapack.MaxColumnSum {
+					work = nanSlice(n)
+				}
+
+				got := impl.Dlanhs(norm, a.Rows, a.Data, lda, work)
+
+				// Zero out A below the first subdiagonal.
+				for i := 2; i < n; i++ {
+					for j := 0; j < max(0, i-1); j++ {
+						a.Data[i*a.Stride+j] = 0
+					}
+				}
+				want := dlange(norm, a.Rows, a.Cols, a.Data, a.Stride)
+
+				if math.Abs(want-got) > tol*want {
+					t.Errorf("Case n=%v,lda=%v,norm=%v: unexpected result. Want %v, got %v.", n, lda, normToString(norm), want, got)
+				}
+			}
+		}
+	}
+}
diff --git a/lapack/testlapack/general.go b/lapack/testlapack/general.go
index 9dd14143..2291dcab 100644
--- a/lapack/testlapack/general.go
+++ b/lapack/testlapack/general.go
@@ -1201,7 +1201,20 @@ func isUpperTriangular(a blas64.General) bool {
 	n := a.Rows
 	for i := 1; i < n; i++ {
 		for j := 0; j < i; j++ {
-			if a.Data[i*a.Stride+j] != 0 {
+			v := a.Data[i*a.Stride+j]
+			if v != 0 || math.IsNaN(v) {
+				return false
+			}
+		}
+	}
+	return true
+}
+
+// isNaNFree returns whether a does not contain NaN elements in reachable elements.
+func isNaNFree(a blas64.General) bool {
+	for i := 0; i < a.Rows; i++ {
+		for j := 0; j < a.Cols; j++ {
+			if math.IsNaN(a.Data[i*a.Stride+j]) {
 				return false
 			}
 		}