lapack/testlapack: clean up hasOrthonormalColumns helper

lapack/testlapack/dorg2l.go

@@ -5,11 +5,9 @@
 package testlapack
 
 import (
-	"math"
 	"testing"
 
 	"golang.org/x/exp/rand"
-
 	"gonum.org/v1/gonum/blas/blas64"
 )
 
@@ -47,30 +45,8 @@ func Dorg2lTest(t *testing.T, impl Dorg2ler) {
 		aCopy := make([]float64, len(a))
 		copy(aCopy, a)
 		impl.Dorg2l(m, n, k, a, lda, tau[n-k:], work)
-		if !hasOrthonormalColumns(m, n, a, lda) {
+		if !hasOrthonormalColumns(blas64.General{m, n, lda, a}) {
 			t.Errorf("Case m=%v, n=%v, k=%v: columns of Q not orthonormal", m, n, k)
 		}
 	}
 }
-
-// hasOrthonormalColumns returns whether the columns of A are orthonormal.
-func hasOrthonormalColumns(m, n int, a []float64, lda int) bool {
-	for i := 0; i < n; i++ {
-		for j := i; j < n; j++ {
-			dot := blas64.Dot(m,
-				blas64.Vector{Inc: lda, Data: a[i:]},
-				blas64.Vector{Inc: lda, Data: a[j:]},
-			)
-			if i == j {
-				if math.Abs(dot-1) > 1e-10 {
-					return false
-				}
-			} else {
-				if math.Abs(dot) > 1e-10 {
-					return false
-				}
-			}
-		}
-	}
-	return true
-}

lapack/testlapack/general.go

@@ -851,6 +851,38 @@ func isOrthogonal(q blas64.General) bool {
 	return true
 }
 
+// hasOrthonormalColumns returns whether the columns of Q are orthonormal.
+func hasOrthonormalColumns(q blas64.General) bool {
+	m, n := q.Rows, q.Cols
+	if m < n {
+		panic("m < n")
+	}
+	ldq := q.Stride
+	const tol = 1e-13
+	for i := 0; i < n; i++ {
+		nrm := blas64.Nrm2(m, blas64.Vector{Data: q.Data[i:], Inc: ldq})
+		if math.IsNaN(nrm) {
+			return false
+		}
+		if math.Abs(nrm-1) > tol {
+			return false
+		}
+		for j := i + 1; j < n; j++ {
+			dot := blas64.Dot(m,
+				blas64.Vector{Data: q.Data[i:], Inc: ldq},
+				blas64.Vector{Data: q.Data[j:], Inc: ldq},
+			)
+			if math.IsNaN(dot) {
+				return false
+			}
+			if math.Abs(dot) > tol {
+				return false
+			}
+		}
+	}
+	return true
+}
+
 // copyMatrix copies an m×n matrix src of stride n into an m×n matrix dst of stride ld.
 func copyMatrix(m, n int, dst []float64, ld int, src []float64) {
 	for i := 0; i < m; i++ {