gonum.org/v1/gonum@v0.14.0/lapack/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  	"fmt"
     9  	"math"
    10  	"testing"
    11  )
    12  
    13  type Dlartger interface {
    14  	Dlartg(f, g float64) (cs, sn, r float64)
    15  }
    16  
    17  func DlartgTest(t *testing.T, impl Dlartger) {
    18  	const tol = 20 * ulp
    19  
    20  	values := []float64{
    21  		-safmax,
    22  		-1 / ulp,
    23  		-1,
    24  		-1.0 / 3,
    25  		-ulp,
    26  		-safmin,
    27  		0,
    28  		safmin,
    29  		ulp,
    30  		1.0 / 3,
    31  		1,
    32  		1 / ulp,
    33  		safmax,
    34  		math.Inf(-1),
    35  		math.Inf(1),
    36  		math.NaN(),
    37  	}
    38  
    39  	for _, f := range values {
    40  		for _, g := range values {
    41  			name := fmt.Sprintf("Case f=%v,g=%v", f, g)
    42  
    43  			// Generate a plane rotation so that
    44  			//  [ cs sn] * [f] = [r]
    45  			//  [-sn cs]   [g] = [0]
    46  			// where cs*cs + sn*sn = 1.
    47  			cs, sn, r := impl.Dlartg(f, g)
    48  
    49  			switch {
    50  			case math.IsNaN(f) || math.IsNaN(g):
    51  				if !math.IsNaN(r) {
    52  					t.Errorf("%v: unexpected r=%v; want NaN", name, r)
    53  				}
    54  			case math.IsInf(f, 0) || math.IsInf(g, 0):
    55  				if !math.IsNaN(r) && !math.IsInf(r, 0) {
    56  					t.Errorf("%v: unexpected r=%v; want NaN or Inf", name, r)
    57  				}
    58  			default:
    59  				d := math.Max(math.Abs(f), math.Abs(g))
    60  				d = math.Min(math.Max(safmin, d), safmax)
    61  				fs := f / d
    62  				gs := g / d
    63  				rs := r / d
    64  
    65  				// Check that cs*f + sn*g = r.
    66  				rnorm := math.Abs(rs)
    67  				if rnorm == 0 {
    68  					rnorm = math.Max(math.Abs(fs), math.Abs(gs))
    69  					if rnorm == 0 {
    70  						rnorm = 1
    71  					}
    72  				}
    73  				resid := math.Abs(rs-(cs*fs+sn*gs)) / rnorm
    74  				if resid > tol {
    75  					t.Errorf("%v: cs*f + sn*g != r; resid=%v", name, resid)
    76  				}
    77  
    78  				// Check that -sn*f + cs*g = 0.
    79  				resid = math.Abs(-sn*fs + cs*gs)
    80  				if resid > tol {
    81  					t.Errorf("%v: -sn*f + cs*g != 0; resid=%v", name, resid)
    82  				}
    83  
    84  				// Check that cs*cs + sn*sn = 1.
    85  				resid = math.Abs(1 - (cs*cs + sn*sn))
    86  				if resid > tol {
    87  					t.Errorf("%v: cs*cs + sn*sn != 1; resid=%v", name, resid)
    88  				}
    89  
    90  				// Check that cs is non-negative.
    91  				if math.Abs(f) > math.Abs(g) && cs < 0 {
    92  					t.Errorf("%v: cs is negative; cs=%v", name, cs)
    93  				}
    94  			}
    95  		}
    96  	}
    97  }