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