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.
10 "golang.org/x/exp/rand"
12 "gonum.org/v1/gonum/blas/blas64"
15 // A123 is the non-symmetric singular matrix
19 // It has three distinct real eigenvalues.
22 func (A123) Matrix() blas64.General {
23 return blas64.General{
35 func (A123) Eigenvalues() []complex128 {
36 return []complex128{16.116843969807043, -1.116843969807043, 0}
39 func (A123) LeftEV() blas64.General {
40 return blas64.General{
45 -0.464547273387671, -0.570795531228578, -0.677043789069485,
46 -0.882905959653586, -0.239520420054206, 0.403865119545174,
47 0.408248290463862, -0.816496580927726, 0.408248290463863,
52 func (A123) RightEV() blas64.General {
53 return blas64.General{
58 -0.231970687246286, -0.785830238742067, 0.408248290463864,
59 -0.525322093301234, -0.086751339256628, -0.816496580927726,
60 -0.818673499356181, 0.612327560228810, 0.408248290463863,
65 // AntisymRandom is a anti-symmetric random matrix. All its eigenvalues are
66 // imaginary with one zero if the order is odd.
67 type AntisymRandom struct {
71 func NewAntisymRandom(n int, rnd *rand.Rand) AntisymRandom {
73 for i := 0; i < n; i++ {
74 for j := i + 1; j < n; j++ {
75 r := rnd.NormFloat64()
76 a.Data[i*a.Stride+j] = r
77 a.Data[j*a.Stride+i] = -r
80 return AntisymRandom{a}
83 func (a AntisymRandom) Matrix() blas64.General {
84 return cloneGeneral(a.mat)
87 func (AntisymRandom) Eigenvalues() []complex128 {
91 // Circulant is a generally non-symmetric matrix given by
92 // A[i,j] = 1 + (j-i+n)%n.
93 // For example, for n=5,
99 // It has real and complex eigenvalues, some possibly repeated.
102 func (c Circulant) Matrix() blas64.General {
105 for i := 0; i < n; i++ {
106 for j := 0; j < n; j++ {
107 a.Data[i*a.Stride+j] = float64(1 + (j-i+n)%n)
113 func (c Circulant) Eigenvalues() []complex128 {
116 ev := make([]complex128, n)
117 for k := 0; k < n; k++ {
118 ev[k] = complex(float64(n), 0)
120 for i := n - 1; i > 0; i-- {
121 for k := 0; k < n; k++ {
122 ev[k] = ev[k]*w[k] + complex(float64(i), 0)
128 // Clement is a generally non-symmetric matrix given by
129 // A[i,j] = i+1, if j == i+1,
130 // = n-i, if j == i-1,
132 // For example, for n=5,
138 // It has n distinct real eigenvalues.
141 func (c Clement) Matrix() blas64.General {
144 for i := 0; i < n; i++ {
146 a.Data[i*a.Stride+i+1] = float64(i + 1)
149 a.Data[i*a.Stride+i-1] = float64(n - i)
155 func (c Clement) Eigenvalues() []complex128 {
157 ev := make([]complex128, n)
159 ev[i] = complex(float64(-n+2*i+1), 0)
164 // Creation is a singular non-symmetric matrix given by
165 // A[i,j] = i, if j == i-1,
167 // For example, for n=5,
173 // Zero is its only eigenvalue.
176 func (c Creation) Matrix() blas64.General {
179 for i := 1; i < n; i++ {
180 a.Data[i*a.Stride+i-1] = float64(i)
185 func (c Creation) Eigenvalues() []complex128 {
186 return make([]complex128, int(c))
189 // Diagonal is a diagonal matrix given by
190 // A[i,j] = i+1, if i == j,
192 // For example, for n=5,
198 // It has n real eigenvalues {1,...,n}.
201 func (d Diagonal) Matrix() blas64.General {
204 for i := 0; i < n; i++ {
205 a.Data[i*a.Stride+i] = float64(i)
210 func (d Diagonal) Eigenvalues() []complex128 {
212 ev := make([]complex128, n)
214 ev[i] = complex(float64(i), 0)
219 // Downshift is a non-singular upper Hessenberg matrix given by
220 // A[i,j] = 1, if (i-j+n)%n == 1,
222 // For example, for n=5,
228 // Its eigenvalues are the complex roots of unity.
231 func (d Downshift) Matrix() blas64.General {
235 for i := 1; i < n; i++ {
236 a.Data[i*a.Stride+i-1] = 1
241 func (d Downshift) Eigenvalues() []complex128 {
242 return rootsOfUnity(int(d))
245 // Fibonacci is an upper Hessenberg matrix with 3 distinct real eigenvalues. For
254 func (f Fibonacci) Matrix() blas64.General {
260 for i := 1; i < n; i++ {
261 a.Data[i*a.Stride+i-1] = 1
262 a.Data[i*a.Stride+i] = 1
267 func (f Fibonacci) Eigenvalues() []complex128 {
269 ev := make([]complex128, n)
270 if n == 0 || n == 1 {
273 phi := 0.5 * (1 + math.Sqrt(5))
274 ev[0] = complex(phi, 0)
275 for i := 1; i < n-1; i++ {
278 ev[n-1] = complex(1-phi, 0)
282 // Gear is a singular non-symmetric matrix with real eigenvalues. For example,
291 func (g Gear) Matrix() blas64.General {
297 for i := 0; i < n-1; i++ {
298 a.Data[i*a.Stride+i+1] = 1
300 for i := 1; i < n; i++ {
301 a.Data[i*a.Stride+i-1] = 1
304 a.Data[(n-1)*a.Stride] = -1
308 func (g Gear) Eigenvalues() []complex128 {
310 ev := make([]complex128, n)
311 if n == 0 || n == 1 {
315 ev[0] = complex(0, 1)
316 ev[1] = complex(0, -1)
323 for p := 1; p <= phi; p++ {
324 ev[w] = complex(float64(2*p)*math.Pi/float64(n), 0)
328 for p := 1; p <= phi; p++ {
329 ev[w] = complex(float64(2*p-1)*math.Pi/float64(n), 0)
332 for i, v := range ev {
333 ev[i] = complex(2*math.Cos(real(v)), 0)
338 // Grcar is an upper Hessenberg matrix given by
339 // A[i,j] = -1 if i == j+1,
340 // = 1 if i <= j and j <= i+k,
342 // For example, for n=5 and k=2,
345 // A = [ . -1 1 1 1 ]
348 // The matrix has sensitive eigenvalues but they are not given explicitly.
354 func (g Grcar) Matrix() blas64.General {
357 for k := 0; k <= g.K; k++ {
358 for i := 0; i < n-k; i++ {
359 a.Data[i*a.Stride+i+k] = 1
362 for i := 1; i < n; i++ {
363 a.Data[i*a.Stride+i-1] = -1
368 func (Grcar) Eigenvalues() []complex128 {
372 // Hanowa is a non-symmetric non-singular matrix of even order given by
373 // A[i,j] = alpha if i == j,
374 // = -i-1 if i < n/2 and j == i + n/2,
375 // = i+1-n/2 if i >= n/2 and j == i - n/2,
377 // The matrix has complex eigenvalues.
379 N int // Order of the matrix, must be even.
383 func (h Hanowa) Matrix() blas64.General {
385 panic("lapack: matrix order must be even")
389 for i := 0; i < n; i++ {
390 a.Data[i*a.Stride+i] = h.Alpha
392 for i := 0; i < n/2; i++ {
393 a.Data[i*a.Stride+i+n/2] = float64(-i - 1)
395 for i := n / 2; i < n; i++ {
396 a.Data[i*a.Stride+i-n/2] = float64(i + 1 - n/2)
401 func (h Hanowa) Eigenvalues() []complex128 {
403 panic("lapack: matrix order must be even")
406 ev := make([]complex128, n)
407 for i := 0; i < n/2; i++ {
408 ev[2*i] = complex(h.Alpha, float64(-i-1))
409 ev[2*i+1] = complex(h.Alpha, float64(i+1))
414 // Lesp is a tridiagonal, generally non-symmetric matrix given by
415 // A[i,j] = -2*i-5 if i == j,
416 // = 1/(i+1) if i == j-1,
417 // = j+1 if i == j+1.
418 // For example, for n=5,
421 // A = [ . 1/3 -9 4 . ]
423 // [ . . . 1/5 -13 ].
424 // The matrix has sensitive eigenvalues but they are not given explicitly.
427 func (l Lesp) Matrix() blas64.General {
430 for i := 0; i < n; i++ {
431 a.Data[i*a.Stride+i] = float64(-2*i - 5)
433 for i := 0; i < n-1; i++ {
434 a.Data[i*a.Stride+i+1] = float64(i + 2)
436 for i := 1; i < n; i++ {
437 a.Data[i*a.Stride+i-1] = 1 / float64(i+1)
442 func (Lesp) Eigenvalues() []complex128 {
446 // Rutis is the 4×4 non-symmetric matrix
451 // It has two distinct real eigenvalues and a pair of complex eigenvalues.
454 func (Rutis) Matrix() blas64.General {
455 return blas64.General{
468 func (Rutis) Eigenvalues() []complex128 {
469 return []complex128{12, 1 + 5i, 1 - 5i, 2}
472 // Tris is a tridiagonal matrix given by
473 // A[i,j] = x if i == j-1,
476 // If x*z is negative, the matrix has complex eigenvalues.
482 func (t Tris) Matrix() blas64.General {
485 for i := 1; i < n; i++ {
486 a.Data[i*a.Stride+i-1] = t.X
488 for i := 0; i < n; i++ {
489 a.Data[i*a.Stride+i] = t.Y
491 for i := 0; i < n-1; i++ {
492 a.Data[i*a.Stride+i+1] = t.Z
497 func (t Tris) Eigenvalues() []complex128 {
499 ev := make([]complex128, n)
501 angle := float64(i+1) * math.Pi / float64(n+1)
504 ev[i] = complex(t.Y+2*math.Sqrt(arg)*math.Cos(angle), 0)
506 ev[i] = complex(t.Y, 2*math.Sqrt(-arg)*math.Cos(angle))
512 // Wilk4 is a 4×4 lower triangular matrix with 4 distinct real eigenvalues.
515 func (Wilk4) Matrix() blas64.General {
516 return blas64.General{
521 0.9143e-4, 0.0, 0.0, 0.0,
522 0.8762, 0.7156e-4, 0.0, 0.0,
523 0.7943, 0.8143, 0.9504e-4, 0.0,
524 0.8017, 0.6123, 0.7165, 0.7123e-4,
529 func (Wilk4) Eigenvalues() []complex128 {
531 0.9504e-4, 0.9143e-4, 0.7156e-4, 0.7123e-4,
535 // Wilk12 is a 12×12 lower Hessenberg matrix with 12 distinct real eigenvalues.
538 func (Wilk12) Matrix() blas64.General {
539 return blas64.General{
544 12, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
545 11, 11, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0,
546 10, 10, 10, 9, 0, 0, 0, 0, 0, 0, 0, 0,
547 9, 9, 9, 9, 8, 0, 0, 0, 0, 0, 0, 0,
548 8, 8, 8, 8, 8, 7, 0, 0, 0, 0, 0, 0,
549 7, 7, 7, 7, 7, 7, 6, 0, 0, 0, 0, 0,
550 6, 6, 6, 6, 6, 6, 6, 5, 0, 0, 0, 0,
551 5, 5, 5, 5, 5, 5, 5, 5, 4, 0, 0, 0,
552 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 0, 0,
553 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 0,
554 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1,
555 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
560 func (Wilk12) Eigenvalues() []complex128 {
577 // Wilk20 is a 20×20 lower Hessenberg matrix. If the parameter is 0, the matrix
578 // has 20 distinct real eigenvalues. If the parameter is 1e-10, the matrix has 6
579 // real eigenvalues and 7 pairs of complex eigenvalues.
582 func (w Wilk20) Matrix() blas64.General {
583 a := zeros(20, 20, 20)
584 for i := 0; i < 20; i++ {
585 a.Data[i*a.Stride+i] = float64(i + 1)
587 for i := 0; i < 19; i++ {
588 a.Data[i*a.Stride+i+1] = 20
590 a.Data[19*a.Stride] = float64(w)
594 func (w Wilk20) Eigenvalues() []complex128 {
596 ev := make([]complex128, 20)
598 ev[i] = complex(float64(i+1), 0)
605 // Zero is a matrix with all elements equal to zero.
608 func (z Zero) Matrix() blas64.General {
610 return zeros(n, n, n)
613 func (z Zero) Eigenvalues() []complex128 {
615 return make([]complex128, n)