gonum.org/v1/gonum@v0.14.0/lapack/testlapack/dlaqr1.go (about)

     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  }