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.
12 "golang.org/x/exp/rand"
14 "gonum.org/v1/gonum/blas"
15 "gonum.org/v1/gonum/blas/blas64"
18 type Dormhrer interface {
19 Dormhr(side blas.Side, trans blas.Transpose, m, n, ilo, ihi int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
24 func DormhrTest(t *testing.T, impl Dormhrer) {
25 rnd := rand.New(rand.NewSource(1))
27 for _, side := range []blas.Side{blas.Left, blas.Right} {
28 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
29 for _, m := range []int{1, 2, 3, 4, 5, 8, 9, 10, 23} {
30 for _, n := range []int{1, 2, 3, 4, 5, 8, 9, 10, 23} {
31 for _, extra := range []int{0, 1, 13} {
32 for cas := 0; cas < 10; cas++ {
34 if side == blas.Right {
42 testDormhr(t, impl, side, trans, m, n, ilo, ihi, extra, true, rnd)
43 testDormhr(t, impl, side, trans, m, n, ilo, ihi, extra, false, rnd)
50 for _, side := range []blas.Side{blas.Left, blas.Right} {
51 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
52 testDormhr(t, impl, side, trans, 0, 0, 0, -1, 0, true, rnd)
53 testDormhr(t, impl, side, trans, 0, 0, 0, -1, 0, false, rnd)
58 func testDormhr(t *testing.T, impl Dormhrer, side blas.Side, trans blas.Transpose, m, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) {
71 // Compute the elementary reflectors and tau.
72 a := randomGeneral(nq, nq, nq+extra, rnd)
75 tau = nanSlice(nq - 1)
77 work := nanSlice(max(1, nq)) // Minimum work for Dgehrd.
78 impl.Dgehrd(nq, ilo, ihi, a.Data, a.Stride, tau, work, len(work))
80 // Construct Q from the elementary reflectors in a and from tau.
83 for j := ilo; j < ihi; j++ {
87 Data: make([]float64, nq),
90 for i := j + 2; i < ihi+1; i++ {
91 v.Data[i] = a.Data[i*a.Stride+j]
93 blas64.Ger(-tau[j], v, v, h)
94 copy(qCopy.Data, q.Data)
95 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
98 c := randomGeneral(m, n, n+extra, rnd)
100 // Compute the product of Q and C explicitly.
101 qc := randomGeneral(m, n, n+extra, rnd)
102 if side == blas.Left {
103 blas64.Gemm(trans, blas.NoTrans, 1, q, c, 0, qc)
105 blas64.Gemm(blas.NoTrans, trans, 1, c, q, 0, qc)
108 // Compute the product of Q and C using Dormhr.
111 impl.Dormhr(side, trans, m, n, ilo, ihi, nil, a.Stride, nil, nil, c.Stride, work, -1)
112 work = nanSlice(int(work[0]))
114 work = nanSlice(max(1, nw))
116 impl.Dormhr(side, trans, m, n, ilo, ihi, a.Data, a.Stride, tau, c.Data, c.Stride, work, len(work))
118 // Compare the two answers.
119 prefix := fmt.Sprintf("Case side=%v, trans=%v, m=%v, n=%v, ilo=%v, ihi=%v, extra=%v, optwork=%v",
120 side, trans, m, n, ilo, ihi, extra, optwork)
121 if !generalOutsideAllNaN(c) {
122 t.Errorf("%v: out-of-range write to C\n%v", prefix, c.Data)
124 for i := 0; i < m; i++ {
125 for j := 0; j < n; j++ {
126 cij := c.Data[i*c.Stride+j]
127 qcij := qc.Data[i*qc.Stride+j]
128 if math.Abs(cij-qcij) > tol {
129 t.Errorf("%v: unexpected value of the QC product at [%v,%v]: want %v, got %v", prefix, i, j, qcij, cij)