1 // Copyright ©2016 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
10 "golang.org/x/exp/rand"
12 "gonum.org/v1/gonum/blas"
13 "gonum.org/v1/gonum/blas/blas64"
14 "gonum.org/v1/gonum/floats"
17 type Dgeql2er interface {
18 Dgeql2(m, n int, a []float64, lda int, tau, work []float64)
21 func Dgeql2Test(t *testing.T, impl Dgeql2er) {
22 rnd := rand.New(rand.NewSource(1))
23 // TODO(btracey): Add tests for m < n.
24 for _, test := range []struct {
37 a := make([]float64, m*lda)
39 a[i] = rnd.NormFloat64()
41 tau := nanSlice(min(m, n))
44 aCopy := make([]float64, len(a))
46 impl.Dgeql2(m, n, a, lda, tau, work)
54 Data: make([]float64, m*m),
56 for i := 0; i < m; i++ {
57 q.Data[i*q.Stride+i] = 1
59 for i := 0; i < k; i++ {
60 h := blas64.General{Rows: m, Cols: m, Stride: m, Data: make([]float64, m*m)}
61 for j := 0; j < m; j++ {
62 h.Data[j*h.Stride+j] = 1
64 v := blas64.Vector{Inc: 1, Data: make([]float64, m)}
66 for j := 0; j < m-k+i; j++ {
67 v.Data[j] = a[j*lda+n-k+i]
69 blas64.Ger(-tau[i], v, v, h)
70 qTmp := blas64.General{Rows: q.Rows, Cols: q.Cols, Stride: q.Stride, Data: make([]float64, len(q.Data))}
71 copy(qTmp.Data, q.Data)
72 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, h, qTmp, 0, q)
74 if !isOrthonormal(q) {
75 t.Errorf("Q is not orthonormal")
81 Data: make([]float64, m*n),
84 for i := m - n; i < m; i++ {
85 for j := 0; j <= min(i-(m-n), n-1); j++ {
86 l.Data[i*l.Stride+j] = a[i*lda+j]
92 ans := blas64.General{Rows: m, Cols: n, Stride: lda, Data: make([]float64, len(a))}
95 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, l, 0, ans)
96 if !floats.EqualApprox(ans.Data, aCopy, 1e-10) {
97 t.Errorf("Reconstruction mismatch: m = %v, n = %v", m, n)