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/floats"
17 type Dormqrer interface {
19 Dormqr(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
22 func DormqrTest(t *testing.T, impl Dormqrer) {
23 rnd := rand.New(rand.NewSource(1))
24 for _, side := range []blas.Side{blas.Left, blas.Right} {
25 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
26 for _, test := range []struct {
27 common, adim, cdim, lda, ldc int
35 {100, 200, 300, 0, 0},
36 {100, 300, 200, 0, 0},
37 {200, 100, 300, 0, 0},
38 {200, 300, 100, 0, 0},
39 {300, 100, 200, 0, 0},
40 {300, 200, 100, 0, 0},
41 {100, 200, 300, 400, 500},
42 {100, 300, 200, 400, 500},
43 {200, 100, 300, 400, 500},
44 {200, 300, 100, 400, 500},
45 {300, 100, 200, 400, 500},
46 {300, 200, 100, 400, 500},
47 {100, 200, 300, 500, 400},
48 {100, 300, 200, 500, 400},
49 {200, 100, 300, 500, 400},
50 {200, 300, 100, 500, 400},
51 {300, 100, 200, 500, 400},
52 {300, 200, 100, 500, 400},
54 var ma, na, mc, nc int
55 if side == blas.Left {
66 // Generate a random matrix
71 a := make([]float64, ma*lda)
75 // Compute random C matrix
80 c := make([]float64, mc*ldc)
87 tau := make([]float64, k)
88 work := make([]float64, 1)
89 impl.Dgeqrf(ma, na, a, lda, tau, work, -1)
90 work = make([]float64, int(work[0]))
91 impl.Dgeqrf(ma, na, a, lda, tau, work, len(work))
93 cCopy := make([]float64, len(c))
95 ans := make([]float64, len(c))
98 if side == blas.Left {
99 work = make([]float64, nc)
101 work = make([]float64, mc)
103 impl.Dorm2r(side, trans, mc, nc, k, a, lda, tau, ans, ldc, work)
105 // Make sure Dorm2r and Dormqr match with small work
106 for i := range work {
107 work[i] = rnd.Float64()
110 impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
111 if !floats.EqualApprox(c, ans, 1e-12) {
112 t.Errorf("Dormqr and Dorm2r mismatch for small work")
115 // Try with the optimum amount of work
117 impl.Dormqr(side, trans, mc, nc, k, nil, lda, nil, nil, ldc, work, -1)
118 work = make([]float64, int(work[0]))
119 for i := range work {
120 work[i] = rnd.Float64()
122 impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
123 if !floats.EqualApprox(c, ans, 1e-12) {
124 t.Errorf("Dormqr and Dorm2r mismatch for full work")
126 for i := 0; i < mc; i++ {
127 fmt.Println(cCopy[i*ldc : (i+1)*ldc])
130 for i := 0; i < mc; i++ {
131 fmt.Println(ans[i*ldc : (i+1)*ldc])
134 for i := 0; i < mc; i++ {
135 fmt.Println(c[i*ldc : (i+1)*ldc])
139 // Try with amount of work that is less than
140 // optimal but still long enough to use the
143 if side == blas.Left {
144 work = make([]float64, 3*nc)
146 work = make([]float64, 3*mc)
148 impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
149 if !floats.EqualApprox(c, ans, 1e-12) {
150 t.Errorf("Dormqr and Dorm2r mismatch for medium work")