vendor/gonum.org/v1/gonum/lapack/testlapack/dlarfx.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 )
  16
  17 type Dlarfxer interface {
  18         Dlarfx(side blas.Side, m, n int, v []float64, tau float64, c []float64, ldc int, work []float64)
  19 }
  20
  21 func DlarfxTest(t *testing.T, impl Dlarfxer) {
  22         rnd := rand.New(rand.NewSource(1))
  23         for _, side := range []blas.Side{blas.Right, blas.Left} {
  24                 for m := 1; m < 12; m++ {
  25                         for n := 1; n < 12; n++ {
  26                                 for _, extra := range []int{0, 1, 11} {
  27                                         for cas := 0; cas < 10; cas++ {
  28                                                 testDlarfx(t, impl, side, m, n, extra, rnd)
  29                                         }
  30                                 }
  31                         }
  32                 }
  33         }
  34 }
  35
  36 func testDlarfx(t *testing.T, impl Dlarfxer, side blas.Side, m, n, extra int, rnd *rand.Rand) {
  37         const tol = 1e-13
  38
  39         c := randomGeneral(m, n, n+extra, rnd)
  40         cWant := randomGeneral(m, n, n+extra, rnd)
  41         tau := rnd.NormFloat64()
  42
  43         var (
  44                 v []float64
  45                 h blas64.General
  46         )
  47         if side == blas.Left {
  48                 v = randomSlice(m, rnd)
  49                 h = eye(m, m+extra)
  50         } else {
  51                 v = randomSlice(n, rnd)
  52                 h = eye(n, n+extra)
  53         }
  54         blas64.Ger(-tau, blas64.Vector{Inc: 1, Data: v}, blas64.Vector{Inc: 1, Data: v}, h)
  55         if side == blas.Left {
  56                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, h, c, 0, cWant)
  57         } else {
  58                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, c, h, 0, cWant)
  59         }
  60
  61         var work []float64
  62         if h.Rows > 10 {
  63                 // Allocate work only if H has order > 10.
  64                 if side == blas.Left {
  65                         work = make([]float64, n)
  66                 } else {
  67                         work = make([]float64, m)
  68                 }
  69         }
  70
  71         impl.Dlarfx(side, m, n, v, tau, c.Data, c.Stride, work)
  72
  73         prefix := fmt.Sprintf("Case side=%v, m=%v, n=%v, extra=%v", side, m, n, extra)
  74
  75         // Check any invalid modifications of c.
  76         if !generalOutsideAllNaN(c) {
  77                 t.Errorf("%v: out-of-range write to C\n%v", prefix, c.Data)
  78         }
  79
  80         if !equalApproxGeneral(c, cWant, tol) {
  81                 t.Errorf("%v: unexpected C\n%v", prefix, c.Data)
  82         }
  83 }