vendor/gonum.org/v1/gonum/lapack/testlapack/dlasy2.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         "math"
  10         "testing"
  11
  12         "golang.org/x/exp/rand"
  13
  14         "gonum.org/v1/gonum/blas"
  15         "gonum.org/v1/gonum/blas/blas64"
  16 )
  17
  18 type Dlasy2er interface {
  19         Dlasy2(tranl, tranr bool, isgn, n1, n2 int, tl []float64, ldtl int, tr []float64, ldtr int, b []float64, ldb int, x []float64, ldx int) (scale, xnorm float64, ok bool)
  20 }
  21
  22 func Dlasy2Test(t *testing.T, impl Dlasy2er) {
  23         rnd := rand.New(rand.NewSource(1))
  24         for _, tranl := range []bool{true, false} {
  25                 for _, tranr := range []bool{true, false} {
  26                         for _, isgn := range []int{1, -1} {
  27                                 for _, n1 := range []int{0, 1, 2} {
  28                                         for _, n2 := range []int{0, 1, 2} {
  29                                                 for _, extra := range []int{0, 1, 2, 13} {
  30                                                         for cas := 0; cas < 1000; cas++ {
  31                                                                 testDlasy2(t, impl, tranl, tranr, isgn, n1, n2, extra, rnd)
  32                                                         }
  33                                                 }
  34                                         }
  35                                 }
  36                         }
  37                 }
  38         }
  39 }
  40
  41 func testDlasy2(t *testing.T, impl Dlasy2er, tranl, tranr bool, isgn, n1, n2, extra int, rnd *rand.Rand) {
  42         const tol = 1e-10
  43
  44         tl := randomGeneral(n1, n1, n1+extra, rnd)
  45         tr := randomGeneral(n2, n2, n2+extra, rnd)
  46         b := randomGeneral(n1, n2, n2+extra, rnd)
  47         x := randomGeneral(n1, n2, n2+extra, rnd)
  48
  49         scale, xnorm, ok := impl.Dlasy2(tranl, tranr, isgn, n1, n2, tl.Data, tl.Stride, tr.Data, tr.Stride, b.Data, b.Stride, x.Data, x.Stride)
  50         if scale > 1 {
  51                 t.Errorf("invalid value of scale, want <= 1, got %v", scale)
  52         }
  53         if n1 == 0 || n2 == 0 {
  54                 return
  55         }
  56
  57         prefix := fmt.Sprintf("Case n1=%v, n2=%v, isgn=%v", n1, n2, isgn)
  58
  59         // Check any invalid modifications of x.
  60         if !generalOutsideAllNaN(x) {
  61                 t.Errorf("%v: out-of-range write to x\n%v", prefix, x.Data)
  62         }
  63
  64         var xnormWant float64
  65         for i := 0; i < n1; i++ {
  66                 var rowsum float64
  67                 for j := 0; j < n2; j++ {
  68                         rowsum += math.Abs(x.Data[i*x.Stride+j])
  69                 }
  70                 if rowsum > xnormWant {
  71                         xnormWant = rowsum
  72                 }
  73         }
  74         if xnormWant != xnorm {
  75                 t.Errorf("%v: unexpected xnorm: want %v, got %v", prefix, xnormWant, xnorm)
  76         }
  77
  78         // Multiply b by scale to get the wanted right-hand side.
  79         for i := 0; i < n1; i++ {
  80                 for j := 0; j < n2; j++ {
  81                         b.Data[i*b.Stride+j] *= scale
  82                 }
  83         }
  84         // Compute the wanted left-hand side.
  85         lhsWant := randomGeneral(n1, n2, n2, rnd)
  86         if tranl {
  87                 blas64.Gemm(blas.Trans, blas.NoTrans, 1, tl, x, 0, lhsWant)
  88         } else {
  89                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tl, x, 0, lhsWant)
  90         }
  91         if tranr {
  92                 blas64.Gemm(blas.NoTrans, blas.Trans, float64(isgn), x, tr, 1, lhsWant)
  93         } else {
  94                 blas64.Gemm(blas.NoTrans, blas.NoTrans, float64(isgn), x, tr, 1, lhsWant)
  95         }
  96         // Compare them.
  97         for i := 0; i < n1; i++ {
  98                 for j := 0; j < n2; j++ {
  99                         diff := lhsWant.Data[i*lhsWant.Stride+j] - b.Data[i*b.Stride+j]
 100                         if math.Abs(diff) > tol && ok {
 101                                 t.Errorf("%v: unexpected result, diff[%v,%v]=%v", prefix, i, j, diff)
 102                         }
 103                 }
 104         }
 105 }