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