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.
17 "golang.org/x/exp/rand"
19 "gonum.org/v1/gonum/blas"
20 "gonum.org/v1/gonum/blas/blas64"
23 type Dlaqr5er interface {
24 Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, nshfts int, sr, si []float64, h []float64, ldh int, iloz, ihiz int, z []float64, ldz int, v []float64, ldv int, u []float64, ldu int, nh int, wh []float64, ldwh int, nv int, wv []float64, ldwv int)
27 type Dlaqr5test struct {
32 ShiftR, ShiftI []float64
39 func Dlaqr5Test(t *testing.T, impl Dlaqr5er) {
40 // Test without using reference data.
41 rnd := rand.New(rand.NewSource(1))
42 for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 30} {
43 for _, extra := range []int{0, 1, 20} {
44 for _, kacc22 := range []int{0, 1, 2} {
45 for cas := 0; cas < 100; cas++ {
46 testDlaqr5(t, impl, n, extra, kacc22, rnd)
52 // Test using reference data computed by the reference netlib
54 file, err := os.Open(filepath.FromSlash("../testlapack/testdata/dlaqr5data.json.gz"))
59 r, err := gzip.NewReader(file)
65 var tests []Dlaqr5test
66 json.NewDecoder(r).Decode(&tests)
67 for _, test := range tests {
70 nshfts := test.NShifts
76 for _, extra := range []int{0, 1, 10} {
77 v := randomGeneral(nshfts/2, 3, 3+extra, rnd)
78 u := randomGeneral(3*nshfts-3, 3*nshfts-3, 3*nshfts-3+extra, rnd)
80 wh := randomGeneral(3*nshfts-3, n, n+extra, rnd)
82 wv := randomGeneral(n, 3*nshfts-3, 3*nshfts-3+extra, rnd)
84 h := nanGeneral(n, n, n+extra)
86 for _, kacc22 := range []int{0, 1, 2} {
87 copyMatrix(n, n, h.Data, h.Stride, test.H)
90 impl.Dlaqr5(wantt, true, kacc22,
94 0, n-1, z.Data, z.Stride,
97 nv, wv.Data, wv.Stride,
98 nh, wh.Data, wh.Stride)
100 prefix := fmt.Sprintf("wantt=%v, n=%v, nshfts=%v, ktop=%v, kbot=%v, extra=%v, kacc22=%v",
101 wantt, n, nshfts, ktop, kbot, extra, kacc22)
102 if !equalApprox(n, n, h.Data, h.Stride, test.HWant, 1e-13) {
103 t.Errorf("Case %v: unexpected matrix H\nh =%v\nhwant=%v", prefix, h.Data, test.HWant)
105 if !equalApprox(n, n, z.Data, z.Stride, test.ZWant, 1e-13) {
106 t.Errorf("Case %v: unexpected matrix Z\nz =%v\nzwant=%v", prefix, z.Data, test.ZWant)
113 func testDlaqr5(t *testing.T, impl Dlaqr5er, n, extra, kacc22 int, rnd *rand.Rand) {
117 sr := make([]float64, nshfts)
118 si := make([]float64, nshfts)
119 for i := 0; i < n; i++ {
120 re := rnd.NormFloat64()
121 im := rnd.NormFloat64()
122 sr[2*i], sr[2*i+1] = re, re
123 si[2*i], si[2*i+1] = im, -im
128 ktop, kbot = kbot, ktop
131 v := randomGeneral(nshfts/2, 3, 3+extra, rnd)
132 u := randomGeneral(3*nshfts-3, 3*nshfts-3, 3*nshfts-3+extra, rnd)
134 wh := randomGeneral(3*nshfts-3, n, n+extra, rnd)
136 wv := randomGeneral(n, 3*nshfts-3, 3*nshfts-3+extra, rnd)
138 h := randomHessenberg(n, n+extra, rnd)
140 h.Data[ktop*h.Stride+ktop-1] = 0
143 h.Data[(kbot+1)*h.Stride+kbot] = 0
146 hCopy.Data = make([]float64, len(h.Data))
147 copy(hCopy.Data, h.Data)
151 impl.Dlaqr5(wantt, wantz, kacc22,
155 0, n-1, z.Data, z.Stride,
158 nv, wv.Data, wv.Stride,
159 nh, wh.Data, wh.Stride)
161 prefix := fmt.Sprintf("Case n=%v, extra=%v, kacc22=%v", n, extra, kacc22)
163 if !generalOutsideAllNaN(h) {
164 t.Errorf("%v: out-of-range write to H\n%v", prefix, h.Data)
166 if !generalOutsideAllNaN(z) {
167 t.Errorf("%v: out-of-range write to Z\n%v", prefix, z.Data)
169 if !generalOutsideAllNaN(u) {
170 t.Errorf("%v: out-of-range write to U\n%v", prefix, u.Data)
172 if !generalOutsideAllNaN(v) {
173 t.Errorf("%v: out-of-range write to V\n%v", prefix, v.Data)
175 if !generalOutsideAllNaN(wh) {
176 t.Errorf("%v: out-of-range write to WH\n%v", prefix, wh.Data)
178 if !generalOutsideAllNaN(wv) {
179 t.Errorf("%v: out-of-range write to WV\n%v", prefix, wv.Data)
182 for i := 0; i < n; i++ {
183 for j := 0; j < i-1; j++ {
184 if h.Data[i*h.Stride+j] != 0 {
185 t.Errorf("%v: H is not Hessenberg, H[%v,%v]!=0", prefix, i, j)
189 if !isOrthonormal(z) {
190 t.Errorf("%v: Z is not orthogonal", prefix)
192 // Construct Z^T * HOrig * Z and check that it is equal to H from Dlaqr5.
193 hz := blas64.General{
197 Data: make([]float64, n*n),
199 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hCopy, z, 0, hz)
200 zhz := blas64.General{
204 Data: make([]float64, n*n),
206 blas64.Gemm(blas.Trans, blas.NoTrans, 1, z, hz, 0, zhz)
207 for i := 0; i < n; i++ {
208 for j := 0; j < n; j++ {
209 diff := zhz.Data[i*zhz.Stride+j] - h.Data[i*h.Stride+j]
210 if math.Abs(diff) > 1e-13 {
211 t.Errorf("%v: Z^T*HOrig*Z and H are not equal, diff at [%v,%v]=%v", prefix, i, j, diff)