test_matrices.go

// Copyright ©2016 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package testlapack

import (
	"math"

	"golang.org/x/exp/rand"

	"gonum.org/v1/gonum/blas/blas64"
)

// A123 is the non-symmetric singular matrix
//
//	    [ 1 2 3 ]
//	A = [ 4 5 6 ]
//	    [ 7 8 9 ]
//
// It has three distinct real eigenvalues.
type A123 struct{}

func (A123) Matrix() blas64.General {
	return blas64.General{
		Rows:   3,
		Cols:   3,
		Stride: 3,
		Data: []float64{
			1, 2, 3,
			4, 5, 6,
			7, 8, 9,
		},
	}
}

func (A123) Eigenvalues() []complex128 {
	return []complex128{16.116843969807043, -1.116843969807043, 0}
}

func (A123) LeftEV() blas64.General {
	return blas64.General{
		Rows:   3,
		Cols:   3,
		Stride: 3,
		Data: []float64{
			-0.464547273387671, -0.570795531228578, -0.677043789069485,
			-0.882905959653586, -0.239520420054206, 0.403865119545174,
			0.408248290463862, -0.816496580927726, 0.408248290463863,
		},
	}
}

func (A123) RightEV() blas64.General {
	return blas64.General{
		Rows:   3,
		Cols:   3,
		Stride: 3,
		Data: []float64{
			-0.231970687246286, -0.785830238742067, 0.408248290463864,
			-0.525322093301234, -0.086751339256628, -0.816496580927726,
			-0.818673499356181, 0.612327560228810, 0.408248290463863,
		},
	}
}

// AntisymRandom is a anti-symmetric random matrix. All its eigenvalues are
// imaginary with one zero if the order is odd.
type AntisymRandom struct {
	mat blas64.General
}

func NewAntisymRandom(n int, rnd *rand.Rand) AntisymRandom {
	a := zeros(n, n, n)
	for i := 0; i < n; i++ {
		for j := i + 1; j < n; j++ {
			r := rnd.NormFloat64()
			a.Data[i*a.Stride+j] = r
			a.Data[j*a.Stride+i] = -r
		}
	}
	return AntisymRandom{a}
}

func (a AntisymRandom) Matrix() blas64.General {
	return cloneGeneral(a.mat)
}

func (AntisymRandom) Eigenvalues() []complex128 {
	return nil
}

// Circulant is a generally non-symmetric matrix given by
//
//	A[i,j] = 1 + (j-i+n)%n.
//
// For example, for n=5,
//
//	    [ 1 2 3 4 5 ]
//	    [ 5 1 2 3 4 ]
//	A = [ 4 5 1 2 3 ]
//	    [ 3 4 5 1 2 ]
//	    [ 2 3 4 5 1 ]
//
// It has real and complex eigenvalues, some possibly repeated.
type Circulant int

func (c Circulant) Matrix() blas64.General {
	n := int(c)
	a := zeros(n, n, n)
	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			a.Data[i*a.Stride+j] = float64(1 + (j-i+n)%n)
		}
	}
	return a
}

func (c Circulant) Eigenvalues() []complex128 {
	n := int(c)
	w := rootsOfUnity(n)
	ev := make([]complex128, n)
	for k := 0; k < n; k++ {
		ev[k] = complex(float64(n), 0)
	}
	for i := n - 1; i > 0; i-- {
		for k := 0; k < n; k++ {
			ev[k] = ev[k]*w[k] + complex(float64(i), 0)
		}
	}
	return ev
}

// Clement is a generally non-symmetric matrix given by
//
//	A[i,j] = i+1  if j == i+1,
//	       = n-i  if j == i-1,
//	       = 0    otherwise.
//
// For example, for n=5,
//
//	    [ . 1 . . . ]
//	    [ 4 . 2 . . ]
//	A = [ . 3 . 3 . ]
//	    [ . . 2 . 4 ]
//	    [ . . . 1 . ]
//
// It has n distinct real eigenvalues.
type Clement int

func (c Clement) Matrix() blas64.General {
	n := int(c)
	a := zeros(n, n, n)
	for i := 0; i < n; i++ {
		if i < n-1 {
			a.Data[i*a.Stride+i+1] = float64(i + 1)
		}
		if i > 0 {
			a.Data[i*a.Stride+i-1] = float64(n - i)
		}
	}
	return a
}

func (c Clement) Eigenvalues() []complex128 {
	n := int(c)
	ev := make([]complex128, n)
	for i := range ev {
		ev[i] = complex(float64(-n+2*i+1), 0)
	}
	return ev
}

// Creation is a singular non-symmetric matrix given by
//
//	A[i,j] = i  if j == i-1,
//	       = 0  otherwise.
//
// For example, for n=5,
//
//	    [ . . . . . ]
//	    [ 1 . . . . ]
//	A = [ . 2 . . . ]
//	    [ . . 3 . . ]
//	    [ . . . 4 . ]
//
// Zero is its only eigenvalue.
type Creation int

func (c Creation) Matrix() blas64.General {
	n := int(c)
	a := zeros(n, n, n)
	for i := 1; i < n; i++ {
		a.Data[i*a.Stride+i-1] = float64(i)
	}
	return a
}

func (c Creation) Eigenvalues() []complex128 {
	return make([]complex128, int(c))
}

// Diagonal is a diagonal matrix given by
//
//	A[i,j] = i+1  if i == j,
//	       = 0    otherwise.
//
// For example, for n=5,
//
//	    [ 1 . . . . ]
//	    [ . 2 . . . ]
//	A = [ . . 3 . . ]
//	    [ . . . 4 . ]
//	    [ . . . . 5 ]
//
// It has n real eigenvalues {1,...,n}.
type Diagonal int

func (d Diagonal) Matrix() blas64.General {
	n := int(d)
	a := zeros(n, n, n)
	for i := 0; i < n; i++ {
		a.Data[i*a.Stride+i] = float64(i)
	}
	return a
}

func (d Diagonal) Eigenvalues() []complex128 {
	n := int(d)
	ev := make([]complex128, n)
	for i := range ev {
		ev[i] = complex(float64(i), 0)
	}
	return ev
}

// Downshift is a non-singular upper Hessenberg matrix given by
//
//	A[i,j] = 1  if (i-j+n)%n == 1,
//	       = 0  otherwise.
//
// For example, for n=5,
//
//	    [ . . . . 1 ]
//	    [ 1 . . . . ]
//	A = [ . 1 . . . ]
//	    [ . . 1 . . ]
//	    [ . . . 1 . ]
//
// Its eigenvalues are the complex roots of unity.
type Downshift int

func (d Downshift) Matrix() blas64.General {
	n := int(d)
	a := zeros(n, n, n)
	a.Data[n-1] = 1
	for i := 1; i < n; i++ {
		a.Data[i*a.Stride+i-1] = 1
	}
	return a
}

func (d Downshift) Eigenvalues() []complex128 {
	return rootsOfUnity(int(d))
}

// Fibonacci is an upper Hessenberg matrix with 3 distinct real eigenvalues. For
// example, for n=5,
//
//	    [ . 1 . . . ]
//	    [ 1 1 . . . ]
//	A = [ . 1 1 . . ]
//	    [ . . 1 1 . ]
//	    [ . . . 1 1 ]
type Fibonacci int

func (f Fibonacci) Matrix() blas64.General {
	n := int(f)
	a := zeros(n, n, n)
	if n > 1 {
		a.Data[1] = 1
	}
	for i := 1; i < n; i++ {
		a.Data[i*a.Stride+i-1] = 1
		a.Data[i*a.Stride+i] = 1
	}
	return a
}

func (f Fibonacci) Eigenvalues() []complex128 {
	n := int(f)
	ev := make([]complex128, n)
	if n == 0 || n == 1 {
		return ev
	}
	phi := 0.5 * (1 + math.Sqrt(5))
	ev[0] = complex(phi, 0)
	for i := 1; i < n-1; i++ {
		ev[i] = 1 + 0i
	}
	ev[n-1] = complex(1-phi, 0)
	return ev
}

// Gear is a singular non-symmetric matrix with real eigenvalues. For example,
// for n=5,
//
//	    [ . 1 . . 1 ]
//	    [ 1 . 1 . . ]
//	A = [ . 1 . 1 . ]
//	    [ . . 1 . 1 ]
//	    [-1 . . 1 . ]
type Gear int

func (g Gear) Matrix() blas64.General {
	n := int(g)
	a := zeros(n, n, n)
	if n == 1 {
		return a
	}
	for i := 0; i < n-1; i++ {
		a.Data[i*a.Stride+i+1] = 1
	}
	for i := 1; i < n; i++ {
		a.Data[i*a.Stride+i-1] = 1
	}
	a.Data[n-1] = 1
	a.Data[(n-1)*a.Stride] = -1
	return a
}

func (g Gear) Eigenvalues() []complex128 {
	n := int(g)
	ev := make([]complex128, n)
	if n == 0 || n == 1 {
		return ev
	}
	if n == 2 {
		ev[0] = complex(0, 1)
		ev[1] = complex(0, -1)
		return ev
	}
	w := 0
	ev[w] = math.Pi / 2
	w++
	phi := (n - 1) / 2
	for p := 1; p <= phi; p++ {
		ev[w] = complex(float64(2*p)*math.Pi/float64(n), 0)
		w++
	}
	phi = n / 2
	for p := 1; p <= phi; p++ {
		ev[w] = complex(float64(2*p-1)*math.Pi/float64(n), 0)
		w++
	}
	for i, v := range ev {
		ev[i] = complex(2*math.Cos(real(v)), 0)
	}
	return ev
}

// Grcar is an upper Hessenberg matrix given by
//
//	A[i,j] = -1  if i == j+1,
//	       = 1   if i <= j and j <= i+k,
//	       = 0   otherwise.
//
// For example, for n=5 and k=2,
//
//	    [  1  1  1  .  . ]
//	    [ -1  1  1  1  . ]
//	A = [  . -1  1  1  1 ]
//	    [  .  . -1  1  1 ]
//	    [  .  .  . -1  1 ]
//
// The matrix has sensitive eigenvalues but they are not given explicitly.
type Grcar struct {
	N int
	K int
}

func (g Grcar) Matrix() blas64.General {
	n := g.N
	a := zeros(n, n, n)
	for k := 0; k <= g.K; k++ {
		for i := 0; i < n-k; i++ {
			a.Data[i*a.Stride+i+k] = 1
		}
	}
	for i := 1; i < n; i++ {
		a.Data[i*a.Stride+i-1] = -1
	}
	return a
}

func (Grcar) Eigenvalues() []complex128 {
	return nil
}

// Hanowa is a non-symmetric non-singular matrix of even order given by
//
//	A[i,j] = alpha    if i == j,
//	       = -i-1     if i < n/2 and j == i + n/2,
//	       = i+1-n/2  if i >= n/2 and j == i - n/2,
//	       = 0        otherwise.
//
// The matrix has complex eigenvalues.
type Hanowa struct {
	N     int // Order of the matrix, must be even.
	Alpha float64
}

func (h Hanowa) Matrix() blas64.General {
	if h.N&0x1 != 0 {
		panic("lapack: matrix order must be even")
	}
	n := h.N
	a := zeros(n, n, n)
	for i := 0; i < n; i++ {
		a.Data[i*a.Stride+i] = h.Alpha
	}
	for i := 0; i < n/2; i++ {
		a.Data[i*a.Stride+i+n/2] = float64(-i - 1)
	}
	for i := n / 2; i < n; i++ {
		a.Data[i*a.Stride+i-n/2] = float64(i + 1 - n/2)
	}
	return a
}

func (h Hanowa) Eigenvalues() []complex128 {
	if h.N&0x1 != 0 {
		panic("lapack: matrix order must be even")
	}
	n := h.N
	ev := make([]complex128, n)
	for i := 0; i < n/2; i++ {
		ev[2*i] = complex(h.Alpha, float64(-i-1))
		ev[2*i+1] = complex(h.Alpha, float64(i+1))
	}
	return ev
}

// Lesp is a tridiagonal, generally non-symmetric matrix given by
//
//	A[i,j] = -2*i-5   if i == j,
//	       = 1/(i+1)  if i == j-1,
//	       = j+1      if i == j+1.
//
// For example, for n=5,
//
//	    [  -5    2    .    .    . ]
//	    [ 1/2   -7    3    .    . ]
//	A = [   .  1/3   -9    4    . ]
//	    [   .    .  1/4  -11    5 ]
//	    [   .    .    .  1/5  -13 ].
//
// The matrix has sensitive eigenvalues but they are not given explicitly.
type Lesp int

func (l Lesp) Matrix() blas64.General {
	n := int(l)
	a := zeros(n, n, n)
	for i := 0; i < n; i++ {
		a.Data[i*a.Stride+i] = float64(-2*i - 5)
	}
	for i := 0; i < n-1; i++ {
		a.Data[i*a.Stride+i+1] = float64(i + 2)
	}
	for i := 1; i < n; i++ {
		a.Data[i*a.Stride+i-1] = 1 / float64(i+1)
	}
	return a
}

func (Lesp) Eigenvalues() []complex128 {
	return nil
}

// Rutis is the 4×4 non-symmetric matrix
//
//	    [ 4 -5  0  3 ]
//	A = [ 0  4 -3 -5 ]
//	    [ 5 -3  4  0 ]
//	    [ 3  0  5  4 ]
//
// It has two distinct real eigenvalues and a pair of complex eigenvalues.
type Rutis struct{}

func (Rutis) Matrix() blas64.General {
	return blas64.General{
		Rows:   4,
		Cols:   4,
		Stride: 4,
		Data: []float64{
			4, -5, 0, 3,
			0, 4, -3, -5,
			5, -3, 4, 0,
			3, 0, 5, 4,
		},
	}
}

func (Rutis) Eigenvalues() []complex128 {
	return []complex128{12, 1 + 5i, 1 - 5i, 2}
}

// Tris is a tridiagonal matrix given by
//
//	A[i,j] = x  if i == j-1,
//	       = y  if i == j,
//	       = z  if i == j+1.
//
// If x*z is negative, the matrix has complex eigenvalues.
type Tris struct {
	N       int
	X, Y, Z float64
}

func (t Tris) Matrix() blas64.General {
	n := t.N
	a := zeros(n, n, n)
	for i := 1; i < n; i++ {
		a.Data[i*a.Stride+i-1] = t.X
	}
	for i := 0; i < n; i++ {
		a.Data[i*a.Stride+i] = t.Y
	}
	for i := 0; i < n-1; i++ {
		a.Data[i*a.Stride+i+1] = t.Z
	}
	return a
}

func (t Tris) Eigenvalues() []complex128 {
	n := t.N
	ev := make([]complex128, n)
	for i := range ev {
		angle := float64(i+1) * math.Pi / float64(n+1)
		arg := t.X * t.Z
		if arg >= 0 {
			ev[i] = complex(t.Y+2*math.Sqrt(arg)*math.Cos(angle), 0)
		} else {
			ev[i] = complex(t.Y, 2*math.Sqrt(-arg)*math.Cos(angle))
		}
	}
	return ev
}

// Wilk4 is a 4×4 lower triangular matrix with 4 distinct real eigenvalues.
type Wilk4 struct{}

func (Wilk4) Matrix() blas64.General {
	return blas64.General{
		Rows:   4,
		Cols:   4,
		Stride: 4,
		Data: []float64{
			0.9143e-4, 0.0, 0.0, 0.0,
			0.8762, 0.7156e-4, 0.0, 0.0,
			0.7943, 0.8143, 0.9504e-4, 0.0,
			0.8017, 0.6123, 0.7165, 0.7123e-4,
		},
	}
}

func (Wilk4) Eigenvalues() []complex128 {
	return []complex128{
		0.9504e-4, 0.9143e-4, 0.7156e-4, 0.7123e-4,
	}
}

// Wilk12 is a 12×12 lower Hessenberg matrix with 12 distinct real eigenvalues.
type Wilk12 struct{}

func (Wilk12) Matrix() blas64.General {
	return blas64.General{
		Rows:   12,
		Cols:   12,
		Stride: 12,
		Data: []float64{
			12, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
			11, 11, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0,
			10, 10, 10, 9, 0, 0, 0, 0, 0, 0, 0, 0,
			9, 9, 9, 9, 8, 0, 0, 0, 0, 0, 0, 0,
			8, 8, 8, 8, 8, 7, 0, 0, 0, 0, 0, 0,
			7, 7, 7, 7, 7, 7, 6, 0, 0, 0, 0, 0,
			6, 6, 6, 6, 6, 6, 6, 5, 0, 0, 0, 0,
			5, 5, 5, 5, 5, 5, 5, 5, 4, 0, 0, 0,
			4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 0, 0,
			3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 0,
			2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1,
			1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
		},
	}
}

func (Wilk12) Eigenvalues() []complex128 {
	return []complex128{
		32.2288915015722210,
		20.1989886458770691,
		12.3110774008685340,
		6.9615330855671154,
		3.5118559485807528,
		1.5539887091319704,
		0.6435053190136506,
		0.2847497205488856,
		0.1436465181918488,
		0.0812276683076552,
		0.0495074140194613,
		0.0310280683208907,
	}
}

// Wilk20 is a 20×20 lower Hessenberg matrix. If the parameter is 0, the matrix
// has 20 distinct real eigenvalues. If the parameter is 1e-10, the matrix has 6
// real eigenvalues and 7 pairs of complex eigenvalues.
type Wilk20 float64

func (w Wilk20) Matrix() blas64.General {
	a := zeros(20, 20, 20)
	for i := 0; i < 20; i++ {
		a.Data[i*a.Stride+i] = float64(i + 1)
	}
	for i := 0; i < 19; i++ {
		a.Data[i*a.Stride+i+1] = 20
	}
	a.Data[19*a.Stride] = float64(w)
	return a
}

func (w Wilk20) Eigenvalues() []complex128 {
	if float64(w) == 0 {
		ev := make([]complex128, 20)
		for i := range ev {
			ev[i] = complex(float64(i+1), 0)
		}
		return ev
	}
	return nil
}

// Zero is a matrix with all elements equal to zero.
type Zero int

func (z Zero) Matrix() blas64.General {
	n := int(z)
	return zeros(n, n, n)
}

func (z Zero) Eigenvalues() []complex128 {
	n := int(z)
	return make([]complex128, n)
}