vendor/gonum.org/v1/gonum/lapack/testlapack/dtrevc3.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"
  10         "testing"
  11
  12         "golang.org/x/exp/rand"
  13
  14         "gonum.org/v1/gonum/blas/blas64"
  15         "gonum.org/v1/gonum/floats"
  16         "gonum.org/v1/gonum/lapack"
  17 )
  18
  19 type Dtrevc3er interface {
  20         Dtrevc3(side lapack.EVSide, howmny lapack.HowMany, selected []bool, n int, t []float64, ldt int, vl []float64, ldvl int, vr []float64, ldvr int, mm int, work []float64, lwork int) int
  21 }
  22
  23 func Dtrevc3Test(t *testing.T, impl Dtrevc3er) {
  24         rnd := rand.New(rand.NewSource(1))
  25         for _, side := range []lapack.EVSide{lapack.RightEV, lapack.LeftEV, lapack.RightLeftEV} {
  26                 for _, howmny := range []lapack.HowMany{lapack.AllEV, lapack.AllEVMulQ, lapack.SelectedEV} {
  27                         for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 34, 100} {
  28                                 for _, extra := range []int{0, 11} {
  29                                         for _, optwork := range []bool{true, false} {
  30                                                 for cas := 0; cas < 10; cas++ {
  31                                                         tmat := randomSchurCanonical(n, n+extra, rnd)
  32                                                         testDtrevc3(t, impl, side, howmny, tmat, optwork, rnd)
  33                                                 }
  34                                         }
  35                                 }
  36                         }
  37                 }
  38         }
  39 }
  40
  41 func testDtrevc3(t *testing.T, impl Dtrevc3er, side lapack.EVSide, howmny lapack.HowMany, tmat blas64.General, optwork bool, rnd *rand.Rand) {
  42         const tol = 1e-14
  43
  44         n := tmat.Rows
  45         extra := tmat.Stride - tmat.Cols
  46         right := side != lapack.LeftEV
  47         left := side != lapack.RightEV
  48
  49         var selected, selectedWant []bool
  50         var mWant int // How many columns will the eigenvectors occupy.
  51         if howmny == lapack.SelectedEV {
  52                 selected = make([]bool, n)
  53                 selectedWant = make([]bool, n)
  54                 // Dtrevc3 will compute only selected eigenvectors. Pick them
  55                 // randomly disregarding whether they are real or complex.
  56                 for i := range selected {
  57                         if rnd.Float64() < 0.5 {
  58                                 selected[i] = true
  59                         }
  60                 }
  61                 // Dtrevc3 will modify (standardize) the slice selected based on
  62                 // whether the corresponding eigenvalues are real or complex. Do
  63                 // the same process here to fill selectedWant.
  64                 for i := 0; i < n; {
  65                         if i == n-1 || tmat.Data[(i+1)*tmat.Stride+i] == 0 {
  66                                 // Real eigenvalue.
  67                                 if selected[i] {
  68                                         selectedWant[i] = true
  69                                         mWant++ // Real eigenvectors occupy one column.
  70                                 }
  71                                 i++
  72                         } else {
  73                                 // Complex eigenvalue.
  74                                 if selected[i] || selected[i+1] {
  75                                         // Dtrevc3 will modify selected so that
  76                                         // only the first element of the pair is
  77                                         // true.
  78                                         selectedWant[i] = true
  79                                         mWant += 2 // Complex eigenvectors occupy two columns.
  80                                 }
  81                                 i += 2
  82                         }
  83                 }
  84         } else {
  85                 // All eigenvectors occupy n columns.
  86                 mWant = n
  87         }
  88
  89         var vr blas64.General
  90         if right {
  91                 if howmny == lapack.AllEVMulQ {
  92                         vr = eye(n, n+extra)
  93                 } else {
  94                         // VR will be overwritten.
  95                         vr = nanGeneral(n, mWant, n+extra)
  96                 }
  97         }
  98
  99         var vl blas64.General
 100         if left {
 101                 if howmny == lapack.AllEVMulQ {
 102                         vl = eye(n, n+extra)
 103                 } else {
 104                         // VL will be overwritten.
 105                         vl = nanGeneral(n, mWant, n+extra)
 106                 }
 107         }
 108
 109         work := make([]float64, max(1, 3*n))
 110         if optwork {
 111                 impl.Dtrevc3(side, howmny, nil, n, nil, 1, nil, 1, nil, 1, mWant, work, -1)
 112                 work = make([]float64, int(work[0]))
 113         }
 114
 115         m := impl.Dtrevc3(side, howmny, selected, n, tmat.Data, tmat.Stride,
 116                 vl.Data, vl.Stride, vr.Data, vr.Stride, mWant, work, len(work))
 117
 118         prefix := fmt.Sprintf("Case side=%v, howmny=%v, n=%v, extra=%v, optwk=%v",
 119                 side, howmny, n, extra, optwork)
 120
 121         if !generalOutsideAllNaN(tmat) {
 122                 t.Errorf("%v: out-of-range write to T", prefix)
 123         }
 124         if !generalOutsideAllNaN(vl) {
 125                 t.Errorf("%v: out-of-range write to VL", prefix)
 126         }
 127         if !generalOutsideAllNaN(vr) {
 128                 t.Errorf("%v: out-of-range write to VR", prefix)
 129         }
 130
 131         if m != mWant {
 132                 t.Errorf("%v: unexpected value of m. Want %v, got %v", prefix, mWant, m)
 133         }
 134
 135         if howmny == lapack.SelectedEV {
 136                 for i := range selected {
 137                         if selected[i] != selectedWant[i] {
 138                                 t.Errorf("%v: unexpected selected[%v]", prefix, i)
 139                         }
 140                 }
 141         }
 142
 143         // Check that the columns of VR and VL are actually eigenvectors and
 144         // that the magnitude of their largest element is 1.
 145         var k int
 146         for j := 0; j < n; {
 147                 re := tmat.Data[j*tmat.Stride+j]
 148                 if j == n-1 || tmat.Data[(j+1)*tmat.Stride+j] == 0 {
 149                         if howmny == lapack.SelectedEV && !selected[j] {
 150                                 j++
 151                                 continue
 152                         }
 153                         if right {
 154                                 ev := columnOf(vr, k)
 155                                 norm := floats.Norm(ev, math.Inf(1))
 156                                 if math.Abs(norm-1) > tol {
 157                                         t.Errorf("%v: magnitude of largest element of VR[:,%v] not 1", prefix, k)
 158                                 }
 159                                 if !isRightEigenvectorOf(tmat, ev, nil, complex(re, 0), tol) {
 160                                         t.Errorf("%v: VR[:,%v] is not real right eigenvector", prefix, k)
 161                                 }
 162                         }
 163                         if left {
 164                                 ev := columnOf(vl, k)
 165                                 norm := floats.Norm(ev, math.Inf(1))
 166                                 if math.Abs(norm-1) > tol {
 167                                         t.Errorf("%v: magnitude of largest element of VL[:,%v] not 1", prefix, k)
 168                                 }
 169                                 if !isLeftEigenvectorOf(tmat, ev, nil, complex(re, 0), tol) {
 170                                         t.Errorf("%v: VL[:,%v] is not real left eigenvector", prefix, k)
 171                                 }
 172                         }
 173                         k++
 174                         j++
 175                         continue
 176                 }
 177                 if howmny == lapack.SelectedEV && !selected[j] {
 178                         j += 2
 179                         continue
 180                 }
 181                 im := math.Sqrt(math.Abs(tmat.Data[(j+1)*tmat.Stride+j])) *
 182                         math.Sqrt(math.Abs(tmat.Data[j*tmat.Stride+j+1]))
 183                 if right {
 184                         evre := columnOf(vr, k)
 185                         evim := columnOf(vr, k+1)
 186                         var evmax float64
 187                         for i, v := range evre {
 188                                 evmax = math.Max(evmax, math.Abs(v)+math.Abs(evim[i]))
 189                         }
 190                         if math.Abs(evmax-1) > tol {
 191                                 t.Errorf("%v: magnitude of largest element of VR[:,%v] not 1", prefix, k)
 192                         }
 193                         if !isRightEigenvectorOf(tmat, evre, evim, complex(re, im), tol) {
 194                                 t.Errorf("%v: VR[:,%v:%v] is not complex right eigenvector", prefix, k, k+1)
 195                         }
 196                         floats.Scale(-1, evim)
 197                         if !isRightEigenvectorOf(tmat, evre, evim, complex(re, -im), tol) {
 198                                 t.Errorf("%v: VR[:,%v:%v] is not complex right eigenvector", prefix, k, k+1)
 199                         }
 200                 }
 201                 if left {
 202                         evre := columnOf(vl, k)
 203                         evim := columnOf(vl, k+1)
 204                         var evmax float64
 205                         for i, v := range evre {
 206                                 evmax = math.Max(evmax, math.Abs(v)+math.Abs(evim[i]))
 207                         }
 208                         if math.Abs(evmax-1) > tol {
 209                                 t.Errorf("%v: magnitude of largest element of VL[:,%v] not 1", prefix, k)
 210                         }
 211                         if !isLeftEigenvectorOf(tmat, evre, evim, complex(re, im), tol) {
 212                                 t.Errorf("%v: VL[:,%v:%v] is not complex left eigenvector", prefix, k, k+1)
 213                         }
 214                         floats.Scale(-1, evim)
 215                         if !isLeftEigenvectorOf(tmat, evre, evim, complex(re, -im), tol) {
 216                                 t.Errorf("%v: VL[:,%v:%v] is not complex left eigenvector", prefix, k, k+1)
 217                         }
 218                 }
 219                 k += 2
 220                 j += 2
 221         }
 222 }