github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/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  	"github.com/jingcheng-WU/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  	const tol = 1e-14
    23  
    24  	rnd := rand.New(rand.NewSource(1))
    25  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
    26  		for _, incx := range []int{1, 3} {
    27  			name := fmt.Sprintf("n=%v,incx=%v", n, incx)
    28  
    29  			// Allocate a slice of minimum length and fill it with
    30  			// random numbers.
    31  			x := make([]float64, max(0, 1+(n-1)*incx))
    32  			for i := range x {
    33  				x[i] = rnd.Float64()
    34  			}
    35  
    36  			// Fill the referenced elements of x and compute the
    37  			// expected result in a non-sophisticated way.
    38  			scale := rnd.Float64()
    39  			ssq := rnd.Float64()
    40  			want := scale * scale * ssq
    41  			for i := 0; i < n; i++ {
    42  				xi := rnd.NormFloat64()
    43  				x[i*incx] = xi
    44  				want += xi * xi
    45  			}
    46  
    47  			xCopy := make([]float64, len(x))
    48  			copy(xCopy, x)
    49  
    50  			// Update scale and ssq so that
    51  			//  scale_out^2 * ssq_out = x[0]^2 + ... + x[n-1]^2 + scale_in^2*ssq_in
    52  			scale, ssq = impl.Dlassq(n, x, incx, scale, ssq)
    53  			if !floats.Equal(x, xCopy) {
    54  				t.Fatalf("%v: unexpected modification of x", name)
    55  			}
    56  
    57  			// Check the result.
    58  			got := scale * scale * ssq
    59  			if math.Abs(got-want) >= tol {
    60  				t.Errorf("%v: unexpected result; got %v, want %v", name, got, want)
    61  			}
    62  		}
    63  	}
    64  }