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.
4
5package testlapack
6
7import (
8	"fmt"
9	"testing"
10
11	"golang.org/x/exp/rand"
12
13	"gonum.org/v1/gonum/blas"
14	"gonum.org/v1/gonum/blas/blas64"
15	"gonum.org/v1/gonum/lapack"
16)
17
18type Dgebaker interface {
19	Dgebak(job lapack.BalanceJob, side lapack.EVSide, n, ilo, ihi int, scale []float64, m int, v []float64, ldv int)
20}
21
22func DgebakTest(t *testing.T, impl Dgebaker) {
23	rnd := rand.New(rand.NewSource(1))
24
25	for _, job := range []lapack.BalanceJob{lapack.BalanceNone, lapack.Permute, lapack.Scale, lapack.PermuteScale} {
26		for _, side := range []lapack.EVSide{lapack.EVLeft, lapack.EVRight} {
27			for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 10, 18, 31, 53} {
28				for _, extra := range []int{0, 11} {
29					for cas := 0; cas < 100; cas++ {
30						m := rnd.Intn(n + 1)
31						v := randomGeneral(n, m, m+extra, rnd)
32						var ilo, ihi int
33						if v.Rows > 0 {
34							ihi = rnd.Intn(n)
35							ilo = rnd.Intn(ihi + 1)
36						} else {
37							ihi = -1
38						}
39						testDgebak(t, impl, job, side, ilo, ihi, v, rnd)
40					}
41				}
42			}
43		}
44	}
45}
46
47func testDgebak(t *testing.T, impl Dgebaker, job lapack.BalanceJob, side lapack.EVSide, ilo, ihi int, v blas64.General, rnd *rand.Rand) {
48	const tol = 1e-15
49	n := v.Rows
50	m := v.Cols
51	extra := v.Stride - v.Cols
52
53	// Create D and D^{-1} by generating random scales between ilo and ihi.
54	d := eye(n, n)
55	dinv := eye(n, n)
56	scale := nanSlice(n)
57	if job == lapack.Scale || job == lapack.PermuteScale {
58		if ilo == ihi {
59			scale[ilo] = 1
60		} else {
61			for i := ilo; i <= ihi; i++ {
62				scale[i] = 2 * rnd.Float64()
63				d.Data[i*d.Stride+i] = scale[i]
64				dinv.Data[i*dinv.Stride+i] = 1 / scale[i]
65			}
66		}
67	}
68
69	// Create P by generating random column swaps.
70	p := eye(n, n)
71	if job == lapack.Permute || job == lapack.PermuteScale {
72		// Make up some random permutations.
73		for i := n - 1; i > ihi; i-- {
74			scale[i] = float64(rnd.Intn(i + 1))
75			blas64.Swap(blas64.Vector{N: n, Data: p.Data[i:], Inc: p.Stride},
76				blas64.Vector{N: n, Data: p.Data[int(scale[i]):], Inc: p.Stride})
77		}
78		for i := 0; i < ilo; i++ {
79			scale[i] = float64(i + rnd.Intn(ihi-i+1))
80			blas64.Swap(blas64.Vector{N: n, Data: p.Data[i:], Inc: p.Stride},
81				blas64.Vector{N: n, Data: p.Data[int(scale[i]):], Inc: p.Stride})
82		}
83	}
84
85	got := cloneGeneral(v)
86	impl.Dgebak(job, side, n, ilo, ihi, scale, m, got.Data, got.Stride)
87
88	prefix := fmt.Sprintf("Case job=%c, side=%c, n=%v, ilo=%v, ihi=%v, m=%v, extra=%v",
89		job, side, n, ilo, ihi, m, extra)
90
91	if !generalOutsideAllNaN(got) {
92		t.Errorf("%v: out-of-range write to V\n%v", prefix, got.Data)
93	}
94
95	// Compute D*V or D^{-1}*V and store into dv.
96	dv := zeros(n, m, m)
97	if side == lapack.EVRight {
98		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, d, v, 0, dv)
99	} else {
100		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, dinv, v, 0, dv)
101	}
102	// Compute P*D*V or P*D^{-1}*V and store into want.
103	want := zeros(n, m, m)
104	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, p, dv, 0, want)
105
106	if !equalApproxGeneral(want, got, tol) {
107		t.Errorf("%v: unexpected value of V", prefix)
108	}
109}
110