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/blas64"
15 "gonum.org/v1/gonum/floats"
16 "gonum.org/v1/gonum/lapack"
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
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)
41 func testDtrevc3(t *testing.T, impl Dtrevc3er, side lapack.EVSide, howmny lapack.HowMany, tmat blas64.General, optwork bool, rnd *rand.Rand) {
45 extra := tmat.Stride - tmat.Cols
46 right := side != lapack.LeftEV
47 left := side != lapack.RightEV
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 {
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.
65 if i == n-1 || tmat.Data[(i+1)*tmat.Stride+i] == 0 {
68 selectedWant[i] = true
69 mWant++ // Real eigenvectors occupy one column.
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
78 selectedWant[i] = true
79 mWant += 2 // Complex eigenvectors occupy two columns.
85 // All eigenvectors occupy n columns.
91 if howmny == lapack.AllEVMulQ {
94 // VR will be overwritten.
95 vr = nanGeneral(n, mWant, n+extra)
101 if howmny == lapack.AllEVMulQ {
104 // VL will be overwritten.
105 vl = nanGeneral(n, mWant, n+extra)
109 work := make([]float64, max(1, 3*n))
111 impl.Dtrevc3(side, howmny, nil, n, nil, 1, nil, 1, nil, 1, mWant, work, -1)
112 work = make([]float64, int(work[0]))
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))
118 prefix := fmt.Sprintf("Case side=%v, howmny=%v, n=%v, extra=%v, optwk=%v",
119 side, howmny, n, extra, optwork)
121 if !generalOutsideAllNaN(tmat) {
122 t.Errorf("%v: out-of-range write to T", prefix)
124 if !generalOutsideAllNaN(vl) {
125 t.Errorf("%v: out-of-range write to VL", prefix)
127 if !generalOutsideAllNaN(vr) {
128 t.Errorf("%v: out-of-range write to VR", prefix)
132 t.Errorf("%v: unexpected value of m. Want %v, got %v", prefix, mWant, m)
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)
143 // Check that the columns of VR and VL are actually eigenvectors and
144 // that the magnitude of their largest element is 1.
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] {
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)
159 if !isRightEigenvectorOf(tmat, ev, nil, complex(re, 0), tol) {
160 t.Errorf("%v: VR[:,%v] is not real right eigenvector", prefix, k)
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)
169 if !isLeftEigenvectorOf(tmat, ev, nil, complex(re, 0), tol) {
170 t.Errorf("%v: VL[:,%v] is not real left eigenvector", prefix, k)
177 if howmny == lapack.SelectedEV && !selected[j] {
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]))
184 evre := columnOf(vr, k)
185 evim := columnOf(vr, k+1)
187 for i, v := range evre {
188 evmax = math.Max(evmax, math.Abs(v)+math.Abs(evim[i]))
190 if math.Abs(evmax-1) > tol {
191 t.Errorf("%v: magnitude of largest element of VR[:,%v] not 1", prefix, k)
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)
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)
202 evre := columnOf(vl, k)
203 evim := columnOf(vl, k+1)
205 for i, v := range evre {
206 evmax = math.Max(evmax, math.Abs(v)+math.Abs(evim[i]))
208 if math.Abs(evmax-1) > tol {
209 t.Errorf("%v: magnitude of largest element of VL[:,%v] not 1", prefix, k)
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)
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)