vendor/gonum.org/v1/gonum/lapack/testlapack/dtrexc.go

   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
   5 package testlapack
   6
   7 import (
   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
  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)
  21 }
  22
  23 func DtrexcTest(t *testing.T, impl Dtrexcer) {
  24         rnd := rand.New(rand.NewSource(1))
  25
  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)
  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.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)
  43                 }
  44         }
  45 }
  46
  47 func testDtrexc(t *testing.T, impl Dtrexcer, compq lapack.EVComp, 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, 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 !isOrthonormal(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 }