github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/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  	"github.com/jingcheng-WU/gonum/blas"
    14  	"github.com/jingcheng-WU/gonum/blas/blas64"
    15  	"github.com/jingcheng-WU/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 < n; i++ {
    44  		re := rnd.NormFloat64()
    45  		im := rnd.NormFloat64()
    46  		sr[2*i], sr[2*i+1] = re, re
    47  		si[2*i], si[2*i+1] = im, -im
    48  	}
    49  	ktop := rnd.Intn(n)
    50  	kbot := rnd.Intn(n)
    51  	if kbot < ktop {
    52  		ktop, kbot = kbot, ktop
    53  	}
    54  
    55  	v := randomGeneral(nshfts/2, 3, 3+extra, rnd)
    56  	u := randomGeneral(3*nshfts-3, 3*nshfts-3, 3*nshfts-3+extra, rnd)
    57  	nh := n
    58  	wh := randomGeneral(3*nshfts-3, n, n+extra, rnd)
    59  	nv := n
    60  	wv := randomGeneral(n, 3*nshfts-3, 3*nshfts-3+extra, rnd)
    61  
    62  	h := randomHessenberg(n, n+extra, rnd)
    63  	if ktop > 0 {
    64  		h.Data[ktop*h.Stride+ktop-1] = 0
    65  	}
    66  	if kbot < n-1 {
    67  		h.Data[(kbot+1)*h.Stride+kbot] = 0
    68  	}
    69  	hCopy := h
    70  	hCopy.Data = make([]float64, len(h.Data))
    71  	copy(hCopy.Data, h.Data)
    72  
    73  	z := eye(n, n+extra)
    74  
    75  	impl.Dlaqr5(wantt, wantz, kacc22,
    76  		n, ktop, kbot,
    77  		nshfts, sr, si,
    78  		h.Data, h.Stride,
    79  		0, n-1, z.Data, z.Stride,
    80  		v.Data, v.Stride,
    81  		u.Data, u.Stride,
    82  		nv, wv.Data, wv.Stride,
    83  		nh, wh.Data, wh.Stride)
    84  
    85  	prefix := fmt.Sprintf("Case n=%v, extra=%v, kacc22=%v", n, extra, kacc22)
    86  
    87  	if !generalOutsideAllNaN(h) {
    88  		t.Errorf("%v: out-of-range write to H\n%v", prefix, h.Data)
    89  	}
    90  	if !generalOutsideAllNaN(z) {
    91  		t.Errorf("%v: out-of-range write to Z\n%v", prefix, z.Data)
    92  	}
    93  	if !generalOutsideAllNaN(u) {
    94  		t.Errorf("%v: out-of-range write to U\n%v", prefix, u.Data)
    95  	}
    96  	if !generalOutsideAllNaN(v) {
    97  		t.Errorf("%v: out-of-range write to V\n%v", prefix, v.Data)
    98  	}
    99  	if !generalOutsideAllNaN(wh) {
   100  		t.Errorf("%v: out-of-range write to WH\n%v", prefix, wh.Data)
   101  	}
   102  	if !generalOutsideAllNaN(wv) {
   103  		t.Errorf("%v: out-of-range write to WV\n%v", prefix, wv.Data)
   104  	}
   105  
   106  	for i := 0; i < n; i++ {
   107  		for j := 0; j < i-1; j++ {
   108  			if h.Data[i*h.Stride+j] != 0 {
   109  				t.Errorf("%v: H is not Hessenberg, H[%v,%v]!=0", prefix, i, j)
   110  			}
   111  		}
   112  	}
   113  	// Check that Z is orthogonal.
   114  	if resid := residualOrthogonal(z, false); resid > tol*float64(n) {
   115  		t.Errorf("Case %v: Z is not orthogonal; resid=%v, want<=%v", prefix, resid, tol*float64(n))
   116  	}
   117  	// Check that |Zᵀ*HOrig*Z - H| is small where H is the result from Dlaqr5.
   118  	hz := zeros(n, n, n)
   119  	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hCopy, z, 0, hz)
   120  	zhz := cloneGeneral(h)
   121  	blas64.Gemm(blas.Trans, blas.NoTrans, 1, z, hz, -1, zhz)
   122  	resid := dlange(lapack.MaxColumnSum, n, n, zhz.Data, zhz.Stride)
   123  	if resid > tol*float64(n) {
   124  		t.Errorf("%v: |Zᵀ*HOrig*Z - H|=%v, want<=%v", prefix, resid, tol*float64(n))
   125  	}
   126  }