vendor/gonum.org/v1/gonum/lapack/testlapack/test_matrices.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         "math"
   9
  10         "golang.org/x/exp/rand"
  11
  12         "gonum.org/v1/gonum/blas/blas64"
  13 )
  14
  15 // A123 is the non-symmetric singular matrix
  16 //      [ 1 2 3 ]
  17 //  A = [ 4 5 6 ]
  18 //      [ 7 8 9 ]
  19 // It has three distinct real eigenvalues.
  20 type A123 struct{}
  21
  22 func (A123) Matrix() blas64.General {
  23         return blas64.General{
  24                 Rows:   3,
  25                 Cols:   3,
  26                 Stride: 3,
  27                 Data: []float64{
  28                         1, 2, 3,
  29                         4, 5, 6,
  30                         7, 8, 9,
  31                 },
  32         }
  33 }
  34
  35 func (A123) Eigenvalues() []complex128 {
  36         return []complex128{16.116843969807043, -1.116843969807043, 0}
  37 }
  38
  39 func (A123) LeftEV() blas64.General {
  40         return blas64.General{
  41                 Rows:   3,
  42                 Cols:   3,
  43                 Stride: 3,
  44                 Data: []float64{
  45                         -0.464547273387671, -0.570795531228578, -0.677043789069485,
  46                         -0.882905959653586, -0.239520420054206, 0.403865119545174,
  47                         0.408248290463862, -0.816496580927726, 0.408248290463863,
  48                 },
  49         }
  50 }
  51
  52 func (A123) RightEV() blas64.General {
  53         return blas64.General{
  54                 Rows:   3,
  55                 Cols:   3,
  56                 Stride: 3,
  57                 Data: []float64{
  58                         -0.231970687246286, -0.785830238742067, 0.408248290463864,
  59                         -0.525322093301234, -0.086751339256628, -0.816496580927726,
  60                         -0.818673499356181, 0.612327560228810, 0.408248290463863,
  61                 },
  62         }
  63 }
  64
  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 {
  68         mat blas64.General
  69 }
  70
  71 func NewAntisymRandom(n int, rnd *rand.Rand) AntisymRandom {
  72         a := zeros(n, n, n)
  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
  78                 }
  79         }
  80         return AntisymRandom{a}
  81 }
  82
  83 func (a AntisymRandom) Matrix() blas64.General {
  84         return cloneGeneral(a.mat)
  85 }
  86
  87 func (AntisymRandom) Eigenvalues() []complex128 {
  88         return nil
  89 }
  90
  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,
  94 //      [ 1 2 3 4 5 ]
  95 //      [ 5 1 2 3 4 ]
  96 //  A = [ 4 5 1 2 3 ]
  97 //      [ 3 4 5 1 2 ]
  98 //      [ 2 3 4 5 1 ]
  99 // It has real and complex eigenvalues, some possibly repeated.
 100 type Circulant int
 101
 102 func (c Circulant) Matrix() blas64.General {
 103         n := int(c)
 104         a := zeros(n, n, n)
 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)
 108                 }
 109         }
 110         return a
 111 }
 112
 113 func (c Circulant) Eigenvalues() []complex128 {
 114         n := int(c)
 115         w := rootsOfUnity(n)
 116         ev := make([]complex128, n)
 117         for k := 0; k < n; k++ {
 118                 ev[k] = complex(float64(n), 0)
 119         }
 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)
 123                 }
 124         }
 125         return ev
 126 }
 127
 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,
 131 //         = 0,    otherwise.
 132 // For example, for n=5,
 133 //      [ . 1 . . . ]
 134 //      [ 4 . 2 . . ]
 135 //  A = [ . 3 . 3 . ]
 136 //      [ . . 2 . 4 ]
 137 //      [ . . . 1 . ]
 138 // It has n distinct real eigenvalues.
 139 type Clement int
 140
 141 func (c Clement) Matrix() blas64.General {
 142         n := int(c)
 143         a := zeros(n, n, n)
 144         for i := 0; i < n; i++ {
 145                 if i < n-1 {
 146                         a.Data[i*a.Stride+i+1] = float64(i + 1)
 147                 }
 148                 if i > 0 {
 149                         a.Data[i*a.Stride+i-1] = float64(n - i)
 150                 }
 151         }
 152         return a
 153 }
 154
 155 func (c Clement) Eigenvalues() []complex128 {
 156         n := int(c)
 157         ev := make([]complex128, n)
 158         for i := range ev {
 159                 ev[i] = complex(float64(-n+2*i+1), 0)
 160         }
 161         return ev
 162 }
 163
 164 // Creation is a singular non-symmetric matrix given by
 165 //  A[i,j] = i,  if j == i-1,
 166 //         = 0,  otherwise.
 167 // For example, for n=5,
 168 //      [ . . . . . ]
 169 //      [ 1 . . . . ]
 170 //  A = [ . 2 . . . ]
 171 //      [ . . 3 . . ]
 172 //      [ . . . 4 . ]
 173 // Zero is its only eigenvalue.
 174 type Creation int
 175
 176 func (c Creation) Matrix() blas64.General {
 177         n := int(c)
 178         a := zeros(n, n, n)
 179         for i := 1; i < n; i++ {
 180                 a.Data[i*a.Stride+i-1] = float64(i)
 181         }
 182         return a
 183 }
 184
 185 func (c Creation) Eigenvalues() []complex128 {
 186         return make([]complex128, int(c))
 187 }
 188
 189 // Diagonal is a diagonal matrix given by
 190 //  A[i,j] = i+1,  if i == j,
 191 //         = 0,    otherwise.
 192 // For example, for n=5,
 193 //      [ 1 . . . . ]
 194 //      [ . 2 . . . ]
 195 //  A = [ . . 3 . . ]
 196 //      [ . . . 4 . ]
 197 //      [ . . . . 5 ]
 198 // It has n real eigenvalues {1,...,n}.
 199 type Diagonal int
 200
 201 func (d Diagonal) Matrix() blas64.General {
 202         n := int(d)
 203         a := zeros(n, n, n)
 204         for i := 0; i < n; i++ {
 205                 a.Data[i*a.Stride+i] = float64(i)
 206         }
 207         return a
 208 }
 209
 210 func (d Diagonal) Eigenvalues() []complex128 {
 211         n := int(d)
 212         ev := make([]complex128, n)
 213         for i := range ev {
 214                 ev[i] = complex(float64(i), 0)
 215         }
 216         return ev
 217 }
 218
 219 // Downshift is a non-singular upper Hessenberg matrix given by
 220 //  A[i,j] = 1,  if (i-j+n)%n == 1,
 221 //         = 0,  otherwise.
 222 // For example, for n=5,
 223 //      [ . . . . 1 ]
 224 //      [ 1 . . . . ]
 225 //  A = [ . 1 . . . ]
 226 //      [ . . 1 . . ]
 227 //      [ . . . 1 . ]
 228 // Its eigenvalues are the complex roots of unity.
 229 type Downshift int
 230
 231 func (d Downshift) Matrix() blas64.General {
 232         n := int(d)
 233         a := zeros(n, n, n)
 234         a.Data[n-1] = 1
 235         for i := 1; i < n; i++ {
 236                 a.Data[i*a.Stride+i-1] = 1
 237         }
 238         return a
 239 }
 240
 241 func (d Downshift) Eigenvalues() []complex128 {
 242         return rootsOfUnity(int(d))
 243 }
 244
 245 // Fibonacci is an upper Hessenberg matrix with 3 distinct real eigenvalues. For
 246 // example, for n=5,
 247 //      [ . 1 . . . ]
 248 //      [ 1 1 . . . ]
 249 //  A = [ . 1 1 . . ]
 250 //      [ . . 1 1 . ]
 251 //      [ . . . 1 1 ]
 252 type Fibonacci int
 253
 254 func (f Fibonacci) Matrix() blas64.General {
 255         n := int(f)
 256         a := zeros(n, n, n)
 257         if n > 1 {
 258                 a.Data[1] = 1
 259         }
 260         for i := 1; i < n; i++ {
 261                 a.Data[i*a.Stride+i-1] = 1
 262                 a.Data[i*a.Stride+i] = 1
 263         }
 264         return a
 265 }
 266
 267 func (f Fibonacci) Eigenvalues() []complex128 {
 268         n := int(f)
 269         ev := make([]complex128, n)
 270         if n == 0 || n == 1 {
 271                 return ev
 272         }
 273         phi := 0.5 * (1 + math.Sqrt(5))
 274         ev[0] = complex(phi, 0)
 275         for i := 1; i < n-1; i++ {
 276                 ev[i] = 1 + 0i
 277         }
 278         ev[n-1] = complex(1-phi, 0)
 279         return ev
 280 }
 281
 282 // Gear is a singular non-symmetric matrix with real eigenvalues. For example,
 283 // for n=5,
 284 //      [ . 1 . . 1 ]
 285 //      [ 1 . 1 . . ]
 286 //  A = [ . 1 . 1 . ]
 287 //      [ . . 1 . 1 ]
 288 //      [-1 . . 1 . ]
 289 type Gear int
 290
 291 func (g Gear) Matrix() blas64.General {
 292         n := int(g)
 293         a := zeros(n, n, n)
 294         if n == 1 {
 295                 return a
 296         }
 297         for i := 0; i < n-1; i++ {
 298                 a.Data[i*a.Stride+i+1] = 1
 299         }
 300         for i := 1; i < n; i++ {
 301                 a.Data[i*a.Stride+i-1] = 1
 302         }
 303         a.Data[n-1] = 1
 304         a.Data[(n-1)*a.Stride] = -1
 305         return a
 306 }
 307
 308 func (g Gear) Eigenvalues() []complex128 {
 309         n := int(g)
 310         ev := make([]complex128, n)
 311         if n == 0 || n == 1 {
 312                 return ev
 313         }
 314         if n == 2 {
 315                 ev[0] = complex(0, 1)
 316                 ev[1] = complex(0, -1)
 317                 return ev
 318         }
 319         w := 0
 320         ev[w] = math.Pi / 2
 321         w++
 322         phi := (n - 1) / 2
 323         for p := 1; p <= phi; p++ {
 324                 ev[w] = complex(float64(2*p)*math.Pi/float64(n), 0)
 325                 w++
 326         }
 327         phi = n / 2
 328         for p := 1; p <= phi; p++ {
 329                 ev[w] = complex(float64(2*p-1)*math.Pi/float64(n), 0)
 330                 w++
 331         }
 332         for i, v := range ev {
 333                 ev[i] = complex(2*math.Cos(real(v)), 0)
 334         }
 335         return ev
 336 }
 337
 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,
 341 //         = 0   otherwise.
 342 // For example, for n=5 and k=2,
 343 //      [  1  1  1  .  . ]
 344 //      [ -1  1  1  1  . ]
 345 //  A = [  . -1  1  1  1 ]
 346 //      [  .  . -1  1  1 ]
 347 //      [  .  .  . -1  1 ]
 348 // The matrix has sensitive eigenvalues but they are not given explicitly.
 349 type Grcar struct {
 350         N int
 351         K int
 352 }
 353
 354 func (g Grcar) Matrix() blas64.General {
 355         n := g.N
 356         a := zeros(n, n, n)
 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
 360                 }
 361         }
 362         for i := 1; i < n; i++ {
 363                 a.Data[i*a.Stride+i-1] = -1
 364         }
 365         return a
 366 }
 367
 368 func (Grcar) Eigenvalues() []complex128 {
 369         return nil
 370 }
 371
 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,
 376 //         = 0        otherwise.
 377 // The matrix has complex eigenvalues.
 378 type Hanowa struct {
 379         N     int // Order of the matrix, must be even.
 380         Alpha float64
 381 }
 382
 383 func (h Hanowa) Matrix() blas64.General {
 384         if h.N&0x1 != 0 {
 385                 panic("lapack: matrix order must be even")
 386         }
 387         n := h.N
 388         a := zeros(n, n, n)
 389         for i := 0; i < n; i++ {
 390                 a.Data[i*a.Stride+i] = h.Alpha
 391         }
 392         for i := 0; i < n/2; i++ {
 393                 a.Data[i*a.Stride+i+n/2] = float64(-i - 1)
 394         }
 395         for i := n / 2; i < n; i++ {
 396                 a.Data[i*a.Stride+i-n/2] = float64(i + 1 - n/2)
 397         }
 398         return a
 399 }
 400
 401 func (h Hanowa) Eigenvalues() []complex128 {
 402         if h.N&0x1 != 0 {
 403                 panic("lapack: matrix order must be even")
 404         }
 405         n := int(h.N)
 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))
 410         }
 411         return ev
 412 }
 413
 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,
 419 //      [  -5    2    .    .    . ]
 420 //      [ 1/2   -7    3    .    . ]
 421 //  A = [   .  1/3   -9    4    . ]
 422 //      [   .    .  1/4  -11    5 ]
 423 //      [   .    .    .  1/5  -13 ].
 424 // The matrix has sensitive eigenvalues but they are not given explicitly.
 425 type Lesp int
 426
 427 func (l Lesp) Matrix() blas64.General {
 428         n := int(l)
 429         a := zeros(n, n, n)
 430         for i := 0; i < n; i++ {
 431                 a.Data[i*a.Stride+i] = float64(-2*i - 5)
 432         }
 433         for i := 0; i < n-1; i++ {
 434                 a.Data[i*a.Stride+i+1] = float64(i + 2)
 435         }
 436         for i := 1; i < n; i++ {
 437                 a.Data[i*a.Stride+i-1] = 1 / float64(i+1)
 438         }
 439         return a
 440 }
 441
 442 func (Lesp) Eigenvalues() []complex128 {
 443         return nil
 444 }
 445
 446 // Rutis is the 4×4 non-symmetric matrix
 447 //      [ 4 -5  0  3 ]
 448 //  A = [ 0  4 -3 -5 ]
 449 //      [ 5 -3  4  0 ]
 450 //      [ 3  0  5  4 ]
 451 // It has two distinct real eigenvalues and a pair of complex eigenvalues.
 452 type Rutis struct{}
 453
 454 func (Rutis) Matrix() blas64.General {
 455         return blas64.General{
 456                 Rows:   4,
 457                 Cols:   4,
 458                 Stride: 4,
 459                 Data: []float64{
 460                         4, -5, 0, 3,
 461                         0, 4, -3, -5,
 462                         5, -3, 4, 0,
 463                         3, 0, 5, 4,
 464                 },
 465         }
 466 }
 467
 468 func (Rutis) Eigenvalues() []complex128 {
 469         return []complex128{12, 1 + 5i, 1 - 5i, 2}
 470 }
 471
 472 // Tris is a tridiagonal matrix given by
 473 //  A[i,j] = x  if i == j-1,
 474 //         = y  if i == j,
 475 //         = z  if i == j+1.
 476 // If x*z is negative, the matrix has complex eigenvalues.
 477 type Tris struct {
 478         N       int
 479         X, Y, Z float64
 480 }
 481
 482 func (t Tris) Matrix() blas64.General {
 483         n := t.N
 484         a := zeros(n, n, n)
 485         for i := 1; i < n; i++ {
 486                 a.Data[i*a.Stride+i-1] = t.X
 487         }
 488         for i := 0; i < n; i++ {
 489                 a.Data[i*a.Stride+i] = t.Y
 490         }
 491         for i := 0; i < n-1; i++ {
 492                 a.Data[i*a.Stride+i+1] = t.Z
 493         }
 494         return a
 495 }
 496
 497 func (t Tris) Eigenvalues() []complex128 {
 498         n := int(t.N)
 499         ev := make([]complex128, n)
 500         for i := range ev {
 501                 angle := float64(i+1) * math.Pi / float64(n+1)
 502                 arg := t.X * t.Z
 503                 if arg >= 0 {
 504                         ev[i] = complex(t.Y+2*math.Sqrt(arg)*math.Cos(angle), 0)
 505                 } else {
 506                         ev[i] = complex(t.Y, 2*math.Sqrt(-arg)*math.Cos(angle))
 507                 }
 508         }
 509         return ev
 510 }
 511
 512 // Wilk4 is a 4×4 lower triangular matrix with 4 distinct real eigenvalues.
 513 type Wilk4 struct{}
 514
 515 func (Wilk4) Matrix() blas64.General {
 516         return blas64.General{
 517                 Rows:   4,
 518                 Cols:   4,
 519                 Stride: 4,
 520                 Data: []float64{
 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,
 525                 },
 526         }
 527 }
 528
 529 func (Wilk4) Eigenvalues() []complex128 {
 530         return []complex128{
 531                 0.9504e-4, 0.9143e-4, 0.7156e-4, 0.7123e-4,
 532         }
 533 }
 534
 535 // Wilk12 is a 12×12 lower Hessenberg matrix with 12 distinct real eigenvalues.
 536 type Wilk12 struct{}
 537
 538 func (Wilk12) Matrix() blas64.General {
 539         return blas64.General{
 540                 Rows:   12,
 541                 Cols:   12,
 542                 Stride: 12,
 543                 Data: []float64{
 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,
 556                 },
 557         }
 558 }
 559
 560 func (Wilk12) Eigenvalues() []complex128 {
 561         return []complex128{
 562                 32.2288915015722210,
 563                 20.1989886458770691,
 564                 12.3110774008685340,
 565                 6.9615330855671154,
 566                 3.5118559485807528,
 567                 1.5539887091319704,
 568                 0.6435053190136506,
 569                 0.2847497205488856,
 570                 0.1436465181918488,
 571                 0.0812276683076552,
 572                 0.0495074140194613,
 573                 0.0310280683208907,
 574         }
 575 }
 576
 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.
 580 type Wilk20 float64
 581
 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)
 586         }
 587         for i := 0; i < 19; i++ {
 588                 a.Data[i*a.Stride+i+1] = 20
 589         }
 590         a.Data[19*a.Stride] = float64(w)
 591         return a
 592 }
 593
 594 func (w Wilk20) Eigenvalues() []complex128 {
 595         if float64(w) == 0 {
 596                 ev := make([]complex128, 20)
 597                 for i := range ev {
 598                         ev[i] = complex(float64(i+1), 0)
 599                 }
 600                 return ev
 601         }
 602         return nil
 603 }
 604
 605 // Zero is a matrix with all elements equal to zero.
 606 type Zero int
 607
 608 func (z Zero) Matrix() blas64.General {
 609         n := int(z)
 610         return zeros(n, n, n)
 611 }
 612
 613 func (z Zero) Eigenvalues() []complex128 {
 614         n := int(z)
 615         return make([]complex128, n)
 616 }