gonum.org/v1/gonum@v0.14.0/lapack/testlapack/dlaqr5.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  	"fmt"
     9  	"testing"
    10  
    11  	"golang.org/x/exp/rand"
    12  
    13  	"gonum.org/v1/gonum/blas"
    14  	"gonum.org/v1/gonum/blas/blas64"
    15  	"gonum.org/v1/gonum/lapack"
    16  )
    17  
    18  type Dlaqr5er interface {
    19  	Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, nshfts int, sr, si []float64, h []float64, ldh int, iloz, ihiz int, z []float64, ldz int, v []float64, ldv int, u []float64, ldu int, nh int, wh []float64, ldwh int, nv int, wv []float64, ldwv int)
    20  }
    21  
    22  func Dlaqr5Test(t *testing.T, impl Dlaqr5er) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 30} {
    25  		for _, extra := range []int{0, 1, 20} {
    26  			for _, kacc22 := range []int{0, 1, 2} {
    27  				for cas := 0; cas < 100; cas++ {
    28  					testDlaqr5(t, impl, n, extra, kacc22, rnd)
    29  				}
    30  			}
    31  		}
    32  	}
    33  }
    34  
    35  func testDlaqr5(t *testing.T, impl Dlaqr5er, n, extra, kacc22 int, rnd *rand.Rand) {
    36  	const tol = 1e-14
    37  
    38  	wantt := true
    39  	wantz := true
    40  	nshfts := 2 * n
    41  	sr := make([]float64, nshfts)
    42  	si := make([]float64, nshfts)
    43  	for i := 0; i < nshfts; {
    44  		if i == nshfts-1 || rnd.Float64() < 0.5 {
    45  			re := rnd.NormFloat64()
    46  			sr[i], si[i] = re, 0
    47  			i++
    48  			continue
    49  		}
    50  		re := rnd.NormFloat64()
    51  		im := rnd.NormFloat64()
    52  		sr[i], sr[i+1] = re, re
    53  		si[i], si[i+1] = im, -im
    54  		i += 2
    55  	}
    56  	ktop := rnd.Intn(n)
    57  	kbot := rnd.Intn(n)
    58  	if kbot < ktop {
    59  		ktop, kbot = kbot, ktop
    60  	}
    61  
    62  	v := randomGeneral(nshfts/2, 3, 3+extra, rnd)
    63  	u := randomGeneral(2*nshfts, 2*nshfts, 2*nshfts+extra, rnd)
    64  	nh := n
    65  	wh := randomGeneral(2*nshfts, n, n+extra, rnd)
    66  	nv := n
    67  	wv := randomGeneral(n, 2*nshfts, 2*nshfts+extra, rnd)
    68  
    69  	h := randomHessenberg(n, n+extra, rnd)
    70  	if ktop > 0 {
    71  		h.Data[ktop*h.Stride+ktop-1] = 0
    72  	}
    73  	if kbot < n-1 {
    74  		h.Data[(kbot+1)*h.Stride+kbot] = 0
    75  	}
    76  	hCopy := h
    77  	hCopy.Data = make([]float64, len(h.Data))
    78  	copy(hCopy.Data, h.Data)
    79  
    80  	z := eye(n, n+extra)
    81  
    82  	impl.Dlaqr5(wantt, wantz, kacc22,
    83  		n, ktop, kbot,
    84  		nshfts, sr, si,
    85  		h.Data, h.Stride,
    86  		0, n-1, z.Data, z.Stride,
    87  		v.Data, v.Stride,
    88  		u.Data, u.Stride,
    89  		nv, wv.Data, wv.Stride,
    90  		nh, wh.Data, wh.Stride)
    91  
    92  	prefix := fmt.Sprintf("Case n=%v, extra=%v, kacc22=%v", n, extra, kacc22)
    93  
    94  	if !generalOutsideAllNaN(h) {
    95  		t.Errorf("%v: out-of-range write to H\n%v", prefix, h.Data)
    96  	}
    97  	if !generalOutsideAllNaN(z) {
    98  		t.Errorf("%v: out-of-range write to Z\n%v", prefix, z.Data)
    99  	}
   100  	if !generalOutsideAllNaN(u) {
   101  		t.Errorf("%v: out-of-range write to U\n%v", prefix, u.Data)
   102  	}
   103  	if !generalOutsideAllNaN(v) {
   104  		t.Errorf("%v: out-of-range write to V\n%v", prefix, v.Data)
   105  	}
   106  	if !generalOutsideAllNaN(wh) {
   107  		t.Errorf("%v: out-of-range write to WH\n%v", prefix, wh.Data)
   108  	}
   109  	if !generalOutsideAllNaN(wv) {
   110  		t.Errorf("%v: out-of-range write to WV\n%v", prefix, wv.Data)
   111  	}
   112  
   113  	for i := 0; i < n; i++ {
   114  		for j := 0; j < i-1; j++ {
   115  			if h.Data[i*h.Stride+j] != 0 {
   116  				t.Errorf("%v: H is not Hessenberg, H[%v,%v]!=0", prefix, i, j)
   117  			}
   118  		}
   119  	}
   120  	// Check that Z is orthogonal.
   121  	if resid := residualOrthogonal(z, false); resid > tol*float64(n) {
   122  		t.Errorf("Case %v: Z is not orthogonal; resid=%v, want<=%v", prefix, resid, tol*float64(n))
   123  	}
   124  	// Check that |Zᵀ*HOrig*Z - H| is small where H is the result from Dlaqr5.
   125  	hz := zeros(n, n, n)
   126  	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hCopy, z, 0, hz)
   127  	zhz := cloneGeneral(h)
   128  	blas64.Gemm(blas.Trans, blas.NoTrans, 1, z, hz, -1, zhz)
   129  	resid := dlange(lapack.MaxColumnSum, n, n, zhz.Data, zhz.Stride)
   130  	if resid > tol*float64(n) {
   131  		t.Errorf("%v: |Zᵀ*HOrig*Z - H|=%v, want<=%v", prefix, resid, tol*float64(n))
   132  	}
   133  }