vendor/gonum.org/v1/gonum/lapack/testlapack/dlaqr1.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         "testing"
  10
  11         "golang.org/x/exp/rand"
  12
  13         "gonum.org/v1/gonum/floats"
  14 )
  15
  16 type Dlaqr1er interface {
  17         Dlaqr1(n int, h []float64, ldh int, sr1, si1, sr2, si2 float64, v []float64)
  18 }
  19
  20 func Dlaqr1Test(t *testing.T, impl Dlaqr1er) {
  21         rnd := rand.New(rand.NewSource(1))
  22
  23         for _, n := range []int{2, 3} {
  24                 for _, ldh := range []int{n, n + 1, n + 10} {
  25                         for _, cas := range []int{1, 2} {
  26                                 for k := 0; k < 100; k++ {
  27                                         v := make([]float64, n)
  28                                         for i := range v {
  29                                                 v[i] = math.NaN()
  30                                         }
  31                                         h := make([]float64, n*(n-1)*ldh)
  32                                         for i := range h {
  33                                                 h[i] = math.NaN()
  34                                         }
  35                                         for i := 0; i < n; i++ {
  36                                                 for j := 0; j < n; j++ {
  37                                                         h[i*ldh+j] = rnd.NormFloat64()
  38                                                 }
  39                                         }
  40                                         var sr1, sr2, si1, si2 float64
  41                                         if cas == 1 {
  42                                                 sr1 = rnd.NormFloat64()
  43                                                 sr2 = sr1
  44                                                 si1 = rnd.NormFloat64()
  45                                                 si2 = -si1
  46                                         } else {
  47                                                 sr1 = rnd.NormFloat64()
  48                                                 sr2 = rnd.NormFloat64()
  49                                                 si1 = 0
  50                                                 si2 = 0
  51                                         }
  52                                         impl.Dlaqr1(n, h, ldh, sr1, si1, sr2, si2, v)
  53
  54                                         // Matrix H - s1*I.
  55                                         h1 := make([]complex128, n*n)
  56                                         for i := 0; i < n; i++ {
  57                                                 for j := 0; j < n; j++ {
  58                                                         h1[i*n+j] = complex(h[i*ldh+j], 0)
  59                                                         if i == j {
  60                                                                 h1[i*n+j] -= complex(sr1, si1)
  61                                                         }
  62                                                 }
  63                                         }
  64                                         // First column of H - s2*I.
  65                                         h2 := make([]complex128, n)
  66                                         for i := 0; i < n; i++ {
  67                                                 h2[i] = complex(h[i*ldh], 0)
  68                                         }
  69                                         h2[0] -= complex(sr2, si2)
  70
  71                                         wantv := make([]float64, n)
  72                                         // Multiply (H-s1*I)*(H-s2*I) to get a tentative
  73                                         // wantv.
  74                                         for i := 0; i < n; i++ {
  75                                                 for j := 0; j < n; j++ {
  76                                                         wantv[i] += real(h1[i*n+j] * h2[j])
  77                                                 }
  78                                         }
  79                                         // Get the unknown scale.
  80                                         scale := v[0] / wantv[0]
  81                                         // Compute the actual wantv.
  82                                         floats.Scale(scale, wantv)
  83
  84                                         // The scale must be the same for all elements.
  85                                         if floats.Distance(wantv, v, math.Inf(1)) > 1e-13 {
  86                                                 t.Errorf("n = %v, ldh = %v, case = %v: Unexpected value of v: got %v, want %v", n, ldh, cas, v, wantv)
  87                                         }
  88                                 }
  89                         }
  90                 }
  91         }
  92 }