package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-10-05 07:06:54 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

testlapack: extend randomSchurCanonical helper

Browse Source 
- generate "bad" matrices with zero and tiny eigenvalues
- return eigenvalues read from the diagonal of generated matrix
This commit is contained in:
Vladimir Chalupecky
2020-02-05 18:41:44 +01:00
committed by Vladimír Chalupecký
parent 5b18cb25e3
commit 1caee46e06
2 changed files with 34 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
lapack/testlapack/general.go
View File

@@ -131,8 +131,8 @@ func randomHessenberg(n, stride int, rnd *rand.Rand) blas64.General {
// form, that is, block upper triangular with 1×1 and 2×2 diagonal blocks where
// each 2×2 diagonal block has its diagonal elements equal and its off-diagonal
// elements of opposite sign.
func randomSchurCanonical(n, stride int, rnd *rand.Rand) blas64.General {
t := randomGeneral(n, n, stride, rnd)
func randomSchurCanonical(n, stride int, rnd *rand.Rand) (t blas64.General, wr, wi []float64) {
t = randomGeneral(n, n, stride, rnd)
// Zero out the lower triangle.
for i := 0; i < t.Rows; i++ {
for j := 0; j < i; j++ {
@@ -141,23 +141,49 @@ func randomSchurCanonical(n, stride int, rnd *rand.Rand) blas64.General {
}
// Randomly create 2×2 diagonal blocks.
for i := 0; i < t.Rows; {
if i == t.Rows-1 || rnd.Float64() < 0.5 {
p := rnd.Float64()
var r float64
switch {
case p < 0.5:
// A half of real parts of eigenvalues will be "normal".
r = t.Data[i*t.Stride+i]
case p < 0.75:
// A quarter will be very small.
r = dlamchS
default:
// A quarter will be zero.
}
// Store the actual value back at the diagonal of T.
t.Data[i*t.Stride+i] = r
// A half of eigenvalues will be real.
if rnd.Float64() < 0.5 || i == t.Rows-1 {
// 1×1 block.
wr = append(wr, r)
wi = append(wi, 0)
i++
continue
}
// 2×2 block.
// Diagonal elements equal.
t.Data[(i+1)*t.Stride+i+1] = t.Data[i*t.Stride+i]
t.Data[(i+1)*t.Stride+i+1] = r
// Element under diagonal very small or "normal".
c := dlamchS
if rnd.Float64() < 0.5 {
c = rnd.NormFloat64()
}
// Off-diagonal elements of opposite sign.
c := rnd.NormFloat64()
if math.Signbit(c) == math.Signbit(t.Data[i*t.Stride+i+1]) {
c *= -1
}
t.Data[(i+1)*t.Stride+i] = c
wr = append(wr, r, r)
c = math.Sqrt(math.Abs(t.Data[i*t.Stride+i+1])) * math.Sqrt(math.Abs(c))
wi = append(wi, c, -c)
i += 2
}
return t
return t, wr, wi
}
// blockedUpperTriGeneral returns a normal random, general matrix in the form