vendor/gonum.org/v1/gonum/lapack/testlapack/dgebak.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         "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 Dgebaker interface {
  19         Dgebak(job lapack.Job, 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.Job{lapack.None, lapack.Permute, lapack.Scale, lapack.PermuteScale} {
  26                 for _, side := range []lapack.EVSide{lapack.LeftEV, lapack.RightEV} {
  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.Job, 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(n,
  76                                 blas64.Vector{p.Stride, p.Data[i:]},
  77                                 blas64.Vector{p.Stride, p.Data[int(scale[i]):]})
  78                 }
  79                 for i := 0; i < ilo; i++ {
  80                         scale[i] = float64(i + rnd.Intn(ihi-i+1))
  81                         blas64.Swap(n,
  82                                 blas64.Vector{p.Stride, p.Data[i:]},
  83                                 blas64.Vector{p.Stride, p.Data[int(scale[i]):]})
  84                 }
  85         }
  86
  87         got := cloneGeneral(v)
  88         impl.Dgebak(job, side, n, ilo, ihi, scale, m, got.Data, got.Stride)
  89
  90         prefix := fmt.Sprintf("Case job=%v, side=%v, n=%v, ilo=%v, ihi=%v, m=%v, extra=%v",
  91                 job, side, n, ilo, ihi, m, extra)
  92
  93         if !generalOutsideAllNaN(got) {
  94                 t.Errorf("%v: out-of-range write to V\n%v", prefix, got.Data)
  95         }
  96
  97         // Compute D*V or D^{-1}*V and store into dv.
  98         dv := zeros(n, m, m)
  99         if side == lapack.RightEV {
 100                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, d, v, 0, dv)
 101         } else {
 102                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, dinv, v, 0, dv)
 103         }
 104         // Compute P*D*V or P*D^{-1}*V and store into want.
 105         want := zeros(n, m, m)
 106         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, p, dv, 0, want)
 107
 108         if !equalApproxGeneral(want, got, tol) {
 109                 t.Errorf("%v: unexpected value of V", prefix)
 110         }
 111 }