diff --git a/lapack/testlapack/dorg2l.go b/lapack/testlapack/dorg2l.go index 62296b71..17f28654 100644 --- a/lapack/testlapack/dorg2l.go +++ b/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 -} diff --git a/lapack/testlapack/general.go b/lapack/testlapack/general.go index 19419bac..cf2e046b 100644 --- a/lapack/testlapack/general.go +++ b/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++ {