github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/testlapack/dlartg.go (about)

     1  // Copyright ©2015 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  	"math"
     9  	"math/rand"
    10  	"testing"
    11  
    12  	"github.com/gonum/floats"
    13  )
    14  
    15  type Dlartger interface {
    16  	Dlartg(f, g float64) (cs, sn, r float64)
    17  }
    18  
    19  func DlartgTest(t *testing.T, impl Dlartger) {
    20  	const tol = 1e-14
    21  	// safmn2 and safmx2 are copied from native.Dlartg.
    22  	safmn2 := math.Pow(dlamchB, math.Trunc(math.Log(dlamchS/dlamchE)/math.Log(dlamchB)/2))
    23  	safmx2 := 1 / safmn2
    24  	rnd := rand.New(rand.NewSource(1))
    25  	for i := 0; i < 1000; i++ {
    26  		var f float64
    27  		var fHuge bool
    28  		switch rnd.Intn(3) {
    29  		case 0:
    30  			// Huge f.
    31  			fHuge = true
    32  			f = math.Pow(10, 10-20*rnd.Float64()) * safmx2
    33  		case 1:
    34  			// Tiny f.
    35  			f = math.Pow(10, 10-20*rnd.Float64()) * safmn2
    36  		default:
    37  			f = rnd.NormFloat64()
    38  		}
    39  		if rnd.Intn(2) == 0 {
    40  			f *= -1
    41  		}
    42  
    43  		var g float64
    44  		var gHuge bool
    45  		switch rnd.Intn(3) {
    46  		case 0:
    47  			// Huge g.
    48  			gHuge = true
    49  			g = math.Pow(10, 10-20*rnd.Float64()) * safmx2
    50  		case 1:
    51  			// Tiny g.
    52  			g = math.Pow(10, 10-20*rnd.Float64()) * safmn2
    53  		default:
    54  			g = rnd.NormFloat64()
    55  		}
    56  		if rnd.Intn(2) == 0 {
    57  			g *= -1
    58  		}
    59  
    60  		cs, sn, r := impl.Dlartg(f, g)
    61  
    62  		rWant := cs*f + sn*g
    63  		if !floats.EqualWithinAbsOrRel(math.Abs(rWant), math.Abs(r), tol, tol) {
    64  			t.Errorf("Case f=%v,g=%v: unexpected r. Want %v, got %v", f, g, rWant, r)
    65  		}
    66  		oneTest := cs*cs + sn*sn
    67  		if math.Abs(oneTest-1) > tol {
    68  			t.Errorf("Case f=%v,g=%v: expected cs^2+sn^2==1, got %v", f, g, oneTest)
    69  		}
    70  		if !fHuge && !gHuge {
    71  			zeroTest := -sn*f + cs*g
    72  			if math.Abs(zeroTest) > tol {
    73  				t.Errorf("Case f=%v,g=%v: expected zero, got %v", f, g, zeroTest)
    74  			}
    75  		}
    76  		if math.Abs(f) > math.Abs(g) && cs < 0 {
    77  			t.Errorf("Case f=%v,g=%v: unexpected negative cs %v", f, g, cs)
    78  		}
    79  	}
    80  	for i := 0; i < 100; i++ {
    81  		cs, sn, _ := impl.Dlartg(rnd.NormFloat64(), 0)
    82  		if cs != 1 {
    83  			t.Errorf("Unexpected cs for g=0. Want 1, got %v", cs)
    84  		}
    85  		if sn != 0 {
    86  			t.Errorf("Unexpected sn for g=0. Want 0, got %v", sn)
    87  		}
    88  	}
    89  	for i := 0; i < 100; i++ {
    90  		cs, sn, _ := impl.Dlartg(0, rnd.NormFloat64())
    91  		if cs != 0 {
    92  			t.Errorf("Unexpected cs for f=0. Want 0, got %v", cs)
    93  		}
    94  		if sn != 1 {
    95  			t.Errorf("Unexpected sn for f=0. Want 1, got %v", sn)
    96  		}
    97  	}
    98  }