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"
16 "gonum.org/v1/gonum/lapack"
19 type Dtrexcer interface {
20 Dtrexc(compq lapack.EVComp, n int, t []float64, ldt int, q []float64, ldq int, ifst, ilst int, work []float64) (ifstOut, ilstOut int, ok bool)
23 func DtrexcTest(t *testing.T, impl Dtrexcer) {
24 rnd := rand.New(rand.NewSource(1))
26 for _, compq := range []lapack.EVComp{lapack.None, 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)
33 testDtrexc(t, impl, compq, tmat, ifst, ilst, extra, rnd)
39 for _, compq := range []lapack.EVComp{lapack.None, 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)
47 func testDtrexc(t *testing.T, impl Dtrexcer, compq lapack.EVComp, tmat blas64.General, ifst, ilst, extra int, rnd *rand.Rand) {
51 fstSize, fstFirst := schurBlockSize(tmat, ifst)
52 lstSize, lstFirst := schurBlockSize(tmat, ilst)
54 tmatCopy := cloneGeneral(tmat)
57 var q, qCopy blas64.General
58 if compq == lapack.UpdateSchur {
61 qCopy = cloneGeneral(q)
66 ifstGot, ilstGot, ok := impl.Dtrexc(compq, n, tmat.Data, tmat.Stride, q.Data, q.Stride, ifst, ilst, work)
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)
71 if !generalOutsideAllNaN(tmat) {
72 t.Errorf("%v: out-of-range write to T", prefix)
74 if wantq && !generalOutsideAllNaN(q) {
75 t.Errorf("%v: out-of-range write to Q", prefix)
79 t.Logf("%v: Dtrexc returned ok=false", prefix)
82 // Check that the index of the first block was correctly updated (if
88 if ifstWant != ifstGot {
89 t.Errorf("%v: unexpected ifst index. Want %v, got %v ", prefix, ifstWant, ifstGot)
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
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.
105 case fstSize == 2 && lstSize == 1:
107 case fstSize == 1 && lstSize == 2:
111 if ilstWant != ilstGot {
112 t.Errorf("%v: unexpected ilst index. Want %v, got %v", prefix, ilstWant, ilstGot)
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)
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)
140 if !isSchurCanonicalGeneral(tmat) {
141 t.Errorf("%v: T is not in Schur canonical form", prefix)
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 {
153 if modMin <= j && j < modMax && j+1 >= i {
156 diff := tmat.Data[i*tmat.Stride+j] - tmatCopy.Data[i*tmatCopy.Stride+j]
158 t.Errorf("%v: unexpected modification at T[%v,%v]", prefix, i, j)
163 // Check that the block at ifstGot was delivered to ilstGot correctly.
165 // 1×1 blocks are swapped exactly.
166 got := tmat.Data[ilstGot*tmat.Stride+ilstGot]
167 want := tmatCopy.Data[ifstGot*tmatCopy.Stride+ifstGot]
169 t.Errorf("%v: unexpected 1×1 block at T[%v,%v]. Want %v, got %v",
170 prefix, want, got, ilstGot, ilstGot)
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)
178 ev1Got, ev2Got := schurBlockEigenvalues(a, b, c, d)
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)
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)
198 if !isOrthonormal(q) {
199 t.Errorf("%v: Q is not orthogonal", prefix)
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 {
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)
212 // Check that Q^T TOrig Q == T.
214 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmatCopy, q, 0, tq)
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)