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.
4
5package testlapack
6
7import (
8	"testing"
9
10	"golang.org/x/exp/rand"
11
12	"gonum.org/v1/gonum/blas"
13	"gonum.org/v1/gonum/blas/blas64"
14	"gonum.org/v1/gonum/floats"
15)
16
17type Dorm2rer interface {
18	Dgeqrfer
19	Dorm2r(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64)
20}
21
22func Dorm2rTest(t *testing.T, impl Dorm2rer) {
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
28			}{
29				{3, 4, 5, 0, 0},
30				{3, 5, 4, 0, 0},
31				{4, 3, 5, 0, 0},
32				{4, 5, 3, 0, 0},
33				{5, 3, 4, 0, 0},
34				{5, 4, 3, 0, 0},
35				{3, 4, 5, 6, 20},
36				{3, 5, 4, 6, 20},
37				{4, 3, 5, 6, 20},
38				{4, 5, 3, 6, 20},
39				{5, 3, 4, 6, 20},
40				{5, 4, 3, 6, 20},
41				{3, 4, 5, 20, 6},
42				{3, 5, 4, 20, 6},
43				{4, 3, 5, 20, 6},
44				{4, 5, 3, 20, 6},
45				{5, 3, 4, 20, 6},
46				{5, 4, 3, 20, 6},
47			} {
48				var ma, na, mc, nc int
49				if side == blas.Left {
50					ma = test.common
51					na = test.adim
52					mc = test.common
53					nc = test.cdim
54				} else {
55					ma = test.common
56					na = test.adim
57					mc = test.cdim
58					nc = test.common
59				}
60
61				// Generate a random matrix
62				lda := test.lda
63				if lda == 0 {
64					lda = na
65				}
66				a := make([]float64, ma*lda)
67				for i := range a {
68					a[i] = rnd.Float64()
69				}
70				ldc := test.ldc
71				if ldc == 0 {
72					ldc = nc
73				}
74				// Compute random C matrix
75				c := make([]float64, mc*ldc)
76				for i := range c {
77					c[i] = rnd.Float64()
78				}
79
80				// Compute QR
81				k := min(ma, na)
82				tau := make([]float64, k)
83				work := make([]float64, 1)
84				impl.Dgeqrf(ma, na, a, lda, tau, work, -1)
85				work = make([]float64, int(work[0]))
86				impl.Dgeqrf(ma, na, a, lda, tau, work, len(work))
87
88				// Build Q from result
89				q := constructQ("QR", ma, na, a, lda, tau)
90
91				cMat := blas64.General{
92					Rows:   mc,
93					Cols:   nc,
94					Stride: ldc,
95					Data:   make([]float64, len(c)),
96				}
97				copy(cMat.Data, c)
98				cMatCopy := blas64.General{
99					Rows:   cMat.Rows,
100					Cols:   cMat.Cols,
101					Stride: cMat.Stride,
102					Data:   make([]float64, len(cMat.Data)),
103				}
104				copy(cMatCopy.Data, cMat.Data)
105				switch {
106				default:
107					panic("bad test")
108				case side == blas.Left && trans == blas.NoTrans:
109					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
110				case side == blas.Left && trans == blas.Trans:
111					blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
112				case side == blas.Right && trans == blas.NoTrans:
113					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMatCopy, q, 0, cMat)
114				case side == blas.Right && trans == blas.Trans:
115					blas64.Gemm(blas.NoTrans, blas.Trans, 1, cMatCopy, q, 0, cMat)
116				}
117				// Do Dorm2r ard compare
118				if side == blas.Left {
119					work = make([]float64, nc)
120				} else {
121					work = make([]float64, mc)
122				}
123				aCopy := make([]float64, len(a))
124				copy(aCopy, a)
125				tauCopy := make([]float64, len(tau))
126				copy(tauCopy, tau)
127				impl.Dorm2r(side, trans, mc, nc, k, a, lda, tau, c, ldc, work)
128				if !floats.Equal(a, aCopy) {
129					t.Errorf("a changed in call")
130				}
131				if !floats.Equal(tau, tauCopy) {
132					t.Errorf("tau changed in call")
133				}
134				if !floats.EqualApprox(cMat.Data, c, 1e-14) {
135					t.Errorf("Multiplication mismatch.\n Want %v \n got %v.", cMat.Data, c)
136				}
137			}
138		}
139	}
140}
141