gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/testlapack/dlassq.go (about)

     1  // Copyright ©2019 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/floats"
    15  )
    16  
    17  type Dlassqer interface {
    18  	Dlassq(n int, x []float64, incx int, scale, ssq float64) (float64, float64)
    19  }
    20  
    21  func DlassqTest(t *testing.T, impl Dlassqer) {
    22  	values := []float64{
    23  		0,
    24  		0.5 * safmin,
    25  		smlnum,
    26  		ulp,
    27  		1,
    28  		1 / ulp,
    29  		bignum,
    30  		safmax,
    31  		math.Inf(1),
    32  		math.NaN(),
    33  	}
    34  
    35  	rnd := rand.New(rand.NewSource(1))
    36  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 30, 40} {
    37  		for _, incx := range []int{1, 3} {
    38  			for cas := 0; cas < 3; cas++ {
    39  				for _, v0 := range values {
    40  					if v0 > 1 {
    41  						v0 *= 0.5
    42  					}
    43  					for _, v1 := range values {
    44  						if v1 > 1 {
    45  							v1 = 0.5 * v1 / math.Sqrt(float64(n+1))
    46  						}
    47  						dlassqTest(t, impl, rnd, n, incx, cas, v0, v1)
    48  					}
    49  				}
    50  			}
    51  		}
    52  	}
    53  }
    54  
    55  func dlassqTest(t *testing.T, impl Dlassqer, rnd *rand.Rand, n, incx, cas int, v0, v1 float64) {
    56  	const (
    57  		rogue = 1234.5678
    58  		tol   = 1e-15
    59  	)
    60  
    61  	name := fmt.Sprintf("n=%v,incx=%v,cas=%v,v0=%v,v1=%v", n, incx, cas, v0, v1)
    62  
    63  	// Generate n random values in (-1,1].
    64  	work := make([]float64, n)
    65  	for i := range work {
    66  		work[i] = 1 - 2*rnd.Float64()
    67  	}
    68  
    69  	// Compute the sum of squares by an unscaled algorithm.
    70  	var workssq float64
    71  	for _, wi := range work {
    72  		workssq += wi * wi
    73  	}
    74  
    75  	// Set initial scale and ssq corresponding to sqrt(v0).
    76  	var scale, ssq float64
    77  	switch cas {
    78  	case 0:
    79  		scale = 1
    80  		ssq = v0
    81  	case 1:
    82  		scale = math.Sqrt(v0)
    83  		ssq = 1
    84  	case 2:
    85  		if v0 < 1 {
    86  			scale = 1.0 / 3
    87  			ssq = 9 * v0
    88  		} else {
    89  			scale = 3
    90  			ssq = v0 / 9
    91  		}
    92  	default:
    93  		panic("bad cas")
    94  	}
    95  
    96  	// Allocate input slice for Dlassq and fill it with
    97  	//  x[:] = scaling factor * work[:].
    98  	x := make([]float64, max(0, 1+(n-1)*incx))
    99  	for i := range x {
   100  		x[i] = rogue
   101  	}
   102  	for i, wi := range work {
   103  		x[i*incx] = v1 * wi
   104  	}
   105  	xCopy := make([]float64, len(x))
   106  	copy(xCopy, x)
   107  
   108  	scaleGot, ssqGot := impl.Dlassq(n, x, incx, scale, ssq)
   109  	nrmGot := scaleGot * math.Sqrt(ssqGot)
   110  
   111  	if !floats.Same(x, xCopy) {
   112  		t.Fatalf("%v: unexpected modification of x", name)
   113  	}
   114  
   115  	// Compute the expected value of the sum of squares.
   116  	z0 := math.Sqrt(v0)
   117  	var z1n float64
   118  	if n >= 1 {
   119  		z1n = v1 * math.Sqrt(workssq)
   120  	}
   121  	zmin := math.Min(z0, z1n)
   122  	zmax := math.Max(z0, z1n)
   123  	var nrmWant float64
   124  	switch {
   125  	case math.IsNaN(z0) || math.IsNaN(z1n):
   126  		nrmWant = math.NaN()
   127  	case zmin == zmax:
   128  		nrmWant = math.Sqrt2 * zmax
   129  	case zmax == 0:
   130  		nrmWant = 0
   131  	default:
   132  		nrmWant = zmax * math.Sqrt(1+(zmin/zmax)*(zmin/zmax))
   133  	}
   134  
   135  	// Check the result.
   136  	switch {
   137  	case math.IsNaN(nrmGot) || math.IsNaN(nrmWant):
   138  		if !math.IsNaN(nrmGot) {
   139  			t.Errorf("%v: expected NaN; got %v", name, nrmGot)
   140  		}
   141  		if !math.IsNaN(nrmWant) {
   142  			t.Errorf("%v: unexpected NaN; want %v", name, nrmWant)
   143  		}
   144  	case nrmGot == nrmWant:
   145  	case nrmWant == 0:
   146  		if nrmGot > tol {
   147  			t.Errorf("%v: unexpected result; got %v, want 0", name, nrmGot)
   148  		}
   149  	default:
   150  		diff := math.Abs(nrmGot-nrmWant) / nrmWant / math.Max(1, float64(n))
   151  		if math.IsNaN(diff) || diff > tol {
   152  			t.Errorf("%v: unexpected result; got %v, want %v", name, nrmGot, nrmWant)
   153  		}
   154  	}
   155  }