1 // Copyright ©2015 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.
11 "golang.org/x/exp/rand"
13 "gonum.org/v1/gonum/blas"
14 "gonum.org/v1/gonum/blas/blas64"
17 type Dgerqfer interface {
18 Dgerqf(m, n int, a []float64, lda int, tau, work []float64, lwork int)
21 func DgerqfTest(t *testing.T, impl Dgerqfer) {
24 rnd := rand.New(rand.NewSource(1))
25 for c, test := range []struct {
61 a := make([]float64, m*lda)
65 aCopy := make([]float64, len(a))
68 tau := make([]float64, k)
70 tau[i] = rnd.Float64()
73 impl.Dgerqf(m, n, a, lda, tau, work, -1)
74 lwkopt := int(work[0])
75 for _, wk := range []struct {
79 {name: "short", length: m},
80 {name: "medium", length: lwkopt - 1},
81 {name: "long", length: lwkopt},
83 if wk.length < max(1, m) {
87 work = make([]float64, lwork)
89 work[i] = rnd.Float64()
92 impl.Dgerqf(m, n, a, lda, tau, work, lwork)
94 // Test that the RQ factorization has completed successfully. Compute
95 // Q based on the vectors.
96 q := constructQ("RQ", m, n, a, lda, tau)
98 // Check that q is orthonormal
99 for i := 0; i < q.Rows; i++ {
100 nrm := blas64.Nrm2(q.Cols, blas64.Vector{Inc: 1, Data: q.Data[i*q.Stride:]})
101 if math.IsNaN(nrm) || math.Abs(nrm-1) > 1e-14 {
102 t.Errorf("Case %v, q not normal", c)
104 for j := 0; j < i; j++ {
105 dot := blas64.Dot(q.Cols, blas64.Vector{Inc: 1, Data: q.Data[i*q.Stride:]}, blas64.Vector{Inc: 1, Data: q.Data[j*q.Stride:]})
106 if math.IsNaN(dot) || math.Abs(dot) > 1e-14 {
107 t.Errorf("Case %v, q not orthogonal", c)
111 // Check that A = R * Q
116 Data: make([]float64, m*n),
118 for i := 0; i < m; i++ {
120 for j := max(0, i-off); j < n; j++ {
121 r.Data[i*r.Stride+j] = a[i*lda+j]
125 got := blas64.General{
129 Data: make([]float64, m*lda),
131 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, r, q, 0, got)
132 want := blas64.General{
138 if !equalApproxGeneral(got, want, tol) {
139 t.Errorf("Case %d, R*Q != a %s\ngot: %+v\nwant:%+v", c, wk.name, got, want)