github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dgebak.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  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/blas/blas64"
    15  	"github.com/gopherd/gonum/lapack"
    16  )
    17  
    18  type Dgebaker interface {
    19  	Dgebak(job lapack.BalanceJob, side lapack.EVSide, n, ilo, ihi int, scale []float64, m int, v []float64, ldv int)
    20  }
    21  
    22  func DgebakTest(t *testing.T, impl Dgebaker) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  
    25  	for _, job := range []lapack.BalanceJob{lapack.BalanceNone, lapack.Permute, lapack.Scale, lapack.PermuteScale} {
    26  		for _, side := range []lapack.EVSide{lapack.EVLeft, lapack.EVRight} {
    27  			for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 10, 18, 31, 53} {
    28  				for _, extra := range []int{0, 11} {
    29  					for cas := 0; cas < 100; cas++ {
    30  						m := rnd.Intn(n + 1)
    31  						v := randomGeneral(n, m, m+extra, rnd)
    32  						var ilo, ihi int
    33  						if v.Rows > 0 {
    34  							ihi = rnd.Intn(n)
    35  							ilo = rnd.Intn(ihi + 1)
    36  						} else {
    37  							ihi = -1
    38  						}
    39  						testDgebak(t, impl, job, side, ilo, ihi, v, rnd)
    40  					}
    41  				}
    42  			}
    43  		}
    44  	}
    45  }
    46  
    47  func testDgebak(t *testing.T, impl Dgebaker, job lapack.BalanceJob, side lapack.EVSide, ilo, ihi int, v blas64.General, rnd *rand.Rand) {
    48  	const tol = 1e-15
    49  	n := v.Rows
    50  	m := v.Cols
    51  	extra := v.Stride - v.Cols
    52  
    53  	// Create D and D^{-1} by generating random scales between ilo and ihi.
    54  	d := eye(n, n)
    55  	dinv := eye(n, n)
    56  	scale := nanSlice(n)
    57  	if job == lapack.Scale || job == lapack.PermuteScale {
    58  		if ilo == ihi {
    59  			scale[ilo] = 1
    60  		} else {
    61  			for i := ilo; i <= ihi; i++ {
    62  				scale[i] = 2 * rnd.Float64()
    63  				d.Data[i*d.Stride+i] = scale[i]
    64  				dinv.Data[i*dinv.Stride+i] = 1 / scale[i]
    65  			}
    66  		}
    67  	}
    68  
    69  	// Create P by generating random column swaps.
    70  	p := eye(n, n)
    71  	if job == lapack.Permute || job == lapack.PermuteScale {
    72  		// Make up some random permutations.
    73  		for i := n - 1; i > ihi; i-- {
    74  			scale[i] = float64(rnd.Intn(i + 1))
    75  			blas64.Swap(blas64.Vector{N: n, Data: p.Data[i:], Inc: p.Stride},
    76  				blas64.Vector{N: n, Data: p.Data[int(scale[i]):], Inc: p.Stride})
    77  		}
    78  		for i := 0; i < ilo; i++ {
    79  			scale[i] = float64(i + rnd.Intn(ihi-i+1))
    80  			blas64.Swap(blas64.Vector{N: n, Data: p.Data[i:], Inc: p.Stride},
    81  				blas64.Vector{N: n, Data: p.Data[int(scale[i]):], Inc: p.Stride})
    82  		}
    83  	}
    84  
    85  	got := cloneGeneral(v)
    86  	impl.Dgebak(job, side, n, ilo, ihi, scale, m, got.Data, got.Stride)
    87  
    88  	prefix := fmt.Sprintf("Case job=%c, side=%c, n=%v, ilo=%v, ihi=%v, m=%v, extra=%v",
    89  		job, side, n, ilo, ihi, m, extra)
    90  
    91  	if !generalOutsideAllNaN(got) {
    92  		t.Errorf("%v: out-of-range write to V\n%v", prefix, got.Data)
    93  	}
    94  
    95  	// Compute D*V or D^{-1}*V and store into dv.
    96  	dv := zeros(n, m, m)
    97  	if side == lapack.EVRight {
    98  		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, d, v, 0, dv)
    99  	} else {
   100  		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, dinv, v, 0, dv)
   101  	}
   102  	// Compute P*D*V or P*D^{-1}*V and store into want.
   103  	want := zeros(n, m, m)
   104  	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, p, dv, 0, want)
   105  
   106  	if !equalApproxGeneral(want, got, tol) {
   107  		t.Errorf("%v: unexpected value of V", prefix)
   108  	}
   109  }