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	"math/cmplx"
10	"testing"
11
12	"golang.org/x/exp/rand"
13
14	"gonum.org/v1/gonum/blas"
15	"gonum.org/v1/gonum/blas/blas64"
16	"gonum.org/v1/gonum/lapack"
17)
18
19type Dtrexcer interface {
20	Dtrexc(compq lapack.UpdateSchurComp, n int, t []float64, ldt int, q []float64, ldq int, ifst, ilst int, work []float64) (ifstOut, ilstOut int, ok bool)
21}
22
23func DtrexcTest(t *testing.T, impl Dtrexcer) {
24	rnd := rand.New(rand.NewSource(1))
25
26	for _, compq := range []lapack.UpdateSchurComp{lapack.UpdateSchurNone, lapack.UpdateSchur} {
27		for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 18, 31, 53} {
28			for _, extra := range []int{0, 1, 11} {
29				for cas := 0; cas < 100; cas++ {
30					tmat := randomSchurCanonical(n, n+extra, rnd)
31					ifst := rnd.Intn(n)
32					ilst := rnd.Intn(n)
33					testDtrexc(t, impl, compq, tmat, ifst, ilst, extra, rnd)
34				}
35			}
36		}
37	}
38
39	for _, compq := range []lapack.UpdateSchurComp{lapack.UpdateSchurNone, lapack.UpdateSchur} {
40		for _, extra := range []int{0, 1, 11} {
41			tmat := randomSchurCanonical(0, extra, rnd)
42			testDtrexc(t, impl, compq, tmat, 0, 0, extra, rnd)
43		}
44	}
45}
46
47func testDtrexc(t *testing.T, impl Dtrexcer, compq lapack.UpdateSchurComp, tmat blas64.General, ifst, ilst, extra int, rnd *rand.Rand) {
48	const tol = 1e-13
49
50	n := tmat.Rows
51	fstSize, fstFirst := schurBlockSize(tmat, ifst)
52	lstSize, lstFirst := schurBlockSize(tmat, ilst)
53
54	tmatCopy := cloneGeneral(tmat)
55
56	var wantq bool
57	var q, qCopy blas64.General
58	if compq == lapack.UpdateSchur {
59		wantq = true
60		q = eye(n, n+extra)
61		qCopy = cloneGeneral(q)
62	}
63
64	work := nanSlice(n)
65
66	ifstGot, ilstGot, ok := impl.Dtrexc(compq, n, tmat.Data, tmat.Stride, q.Data, max(1, q.Stride), ifst, ilst, work)
67
68	prefix := fmt.Sprintf("Case compq=%v, n=%v, ifst=%v, nbf=%v, ilst=%v, nbl=%v, extra=%v",
69		compq, n, ifst, fstSize, ilst, lstSize, extra)
70
71	if !generalOutsideAllNaN(tmat) {
72		t.Errorf("%v: out-of-range write to T", prefix)
73	}
74	if wantq && !generalOutsideAllNaN(q) {
75		t.Errorf("%v: out-of-range write to Q", prefix)
76	}
77
78	if !ok {
79		t.Logf("%v: Dtrexc returned ok=false", prefix)
80	}
81
82	// Check that the index of the first block was correctly updated (if
83	// necessary).
84	ifstWant := ifst
85	if !fstFirst {
86		ifstWant = ifst - 1
87	}
88	if ifstWant != ifstGot {
89		t.Errorf("%v: unexpected ifst index. Want %v, got %v ", prefix, ifstWant, ifstGot)
90	}
91
92	// Check that the index of the last block is as expected when ok=true.
93	// When ok=false, we don't know at which block the algorithm failed, so
94	// we don't check.
95	ilstWant := ilst
96	if !lstFirst {
97		ilstWant--
98	}
99	if ok {
100		if ifstWant < ilstWant {
101			// If the blocks are swapped backwards, these
102			// adjustments are not necessary, the first row of the
103			// last block will end up at ifst.
104			switch {
105			case fstSize == 2 && lstSize == 1:
106				ilstWant--
107			case fstSize == 1 && lstSize == 2:
108				ilstWant++
109			}
110		}
111		if ilstWant != ilstGot {
112			t.Errorf("%v: unexpected ilst index. Want %v, got %v", prefix, ilstWant, ilstGot)
113		}
114	}
115
116	if n <= 1 || ifstGot == ilstGot {
117		// Too small matrix or no swapping.
118		// Check that T was not modified.
119		for i := 0; i < n; i++ {
120			for j := 0; j < n; j++ {
121				if tmat.Data[i*tmat.Stride+j] != tmatCopy.Data[i*tmatCopy.Stride+j] {
122					t.Errorf("%v: unexpected modification at T[%v,%v]", prefix, i, j)
123				}
124			}
125		}
126		if !wantq {
127			return
128		}
129		// Check that Q was not modified.
130		for i := 0; i < n; i++ {
131			for j := 0; j < n; j++ {
132				if q.Data[i*q.Stride+j] != qCopy.Data[i*qCopy.Stride+j] {
133					t.Errorf("%v: unexpected modification at Q[%v,%v]", prefix, i, j)
134				}
135			}
136		}
137		return
138	}
139
140	if !isSchurCanonicalGeneral(tmat) {
141		t.Errorf("%v: T is not in Schur canonical form", prefix)
142	}
143
144	// Check that T was not modified except above the second subdiagonal in
145	// rows and columns [modMin,modMax].
146	modMin := min(ifstGot, ilstGot)
147	modMax := max(ifstGot, ilstGot) + fstSize
148	for i := 0; i < n; i++ {
149		for j := 0; j < n; j++ {
150			if modMin <= i && i < modMax && j+1 >= i {
151				continue
152			}
153			if modMin <= j && j < modMax && j+1 >= i {
154				continue
155			}
156			diff := tmat.Data[i*tmat.Stride+j] - tmatCopy.Data[i*tmatCopy.Stride+j]
157			if diff != 0 {
158				t.Errorf("%v: unexpected modification at T[%v,%v]", prefix, i, j)
159			}
160		}
161	}
162
163	// Check that the block at ifstGot was delivered to ilstGot correctly.
164	if fstSize == 1 {
165		// 1×1 blocks are swapped exactly.
166		got := tmat.Data[ilstGot*tmat.Stride+ilstGot]
167		want := tmatCopy.Data[ifstGot*tmatCopy.Stride+ifstGot]
168		if want != got {
169			t.Errorf("%v: unexpected 1×1 block at T[%v,%v]. Want %v, got %v",
170				prefix, want, got, ilstGot, ilstGot)
171		}
172	} else {
173		// Check that the swapped 2×2 block is in Schur canonical form.
174		a, b, c, d := extract2x2Block(tmat.Data[ilstGot*tmat.Stride+ilstGot:], tmat.Stride)
175		if !isSchurCanonical(a, b, c, d) {
176			t.Errorf("%v: 2×2 block at T[%v,%v] not in Schur canonical form", prefix, ilstGot, ilstGot)
177		}
178		ev1Got, ev2Got := schurBlockEigenvalues(a, b, c, d)
179
180		// Check that the swapped 2×2 block has the same eigenvalues.
181		// The block was originally located at T[ifstGot,ifstGot].
182		a, b, c, d = extract2x2Block(tmatCopy.Data[ifstGot*tmatCopy.Stride+ifstGot:], tmatCopy.Stride)
183		ev1Want, ev2Want := schurBlockEigenvalues(a, b, c, d)
184		if cmplx.Abs(ev1Got-ev1Want) > tol {
185			t.Errorf("%v: unexpected first eigenvalue of 2×2 block at T[%v,%v]. Want %v, got %v",
186				prefix, ilstGot, ilstGot, ev1Want, ev1Got)
187		}
188		if cmplx.Abs(ev2Got-ev2Want) > tol {
189			t.Errorf("%v: unexpected second eigenvalue of 2×2 block at T[%v,%v]. Want %v, got %v",
190				prefix, ilstGot, ilstGot, ev2Want, ev2Got)
191		}
192	}
193
194	if !wantq {
195		return
196	}
197
198	if !isOrthogonal(q) {
199		t.Errorf("%v: Q is not orthogonal", prefix)
200	}
201	// Check that Q is unchanged outside of columns [modMin,modMax].
202	for i := 0; i < n; i++ {
203		for j := 0; j < n; j++ {
204			if modMin <= j && j < modMax {
205				continue
206			}
207			if q.Data[i*q.Stride+j]-qCopy.Data[i*qCopy.Stride+j] != 0 {
208				t.Errorf("%v: unexpected modification of Q[%v,%v]", prefix, i, j)
209			}
210		}
211	}
212	// Check that Q^T TOrig Q == T.
213	tq := eye(n, n)
214	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmatCopy, q, 0, tq)
215	qtq := eye(n, n)
216	blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, tq, 0, qtq)
217	if !equalApproxGeneral(qtq, tmat, tol) {
218		t.Errorf("%v: Q^T (initial T) Q and (final T) are not equal", prefix)
219	}
220}
221