github.com/gopherd/gonum@v0.0.4/lapack/gonum/dlassq.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 gonum
     6  
     7  import "math"
     8  
     9  // Dlassq updates a sum of squares represented in scaled form. Dlassq returns
    10  // the values scl and smsq such that
    11  //  scl^2*smsq = X[0]^2 + ... + X[n-1]^2 + scale^2*sumsq
    12  // The value of sumsq is assumed to be non-negative.
    13  //
    14  // Dlassq is an internal routine. It is exported for testing purposes.
    15  func (impl Implementation) Dlassq(n int, x []float64, incx int, scale float64, sumsq float64) (scl, smsq float64) {
    16  	// Implementation based on Supplemental Material to:
    17  	// Edward Anderson. 2017. Algorithm 978: Safe Scaling in the Level 1 BLAS.
    18  	// ACM Trans. Math. Softw. 44, 1, Article 12 (July 2017), 28 pages.
    19  	// DOI: https://doi.org/10.1145/3061665
    20  	switch {
    21  	case n < 0:
    22  		panic(nLT0)
    23  	case incx <= 0:
    24  		panic(badIncX)
    25  	case len(x) < 1+(n-1)*incx:
    26  		panic(shortX)
    27  	}
    28  
    29  	if math.IsNaN(scale) || math.IsNaN(sumsq) {
    30  		return scale, sumsq
    31  	}
    32  
    33  	if sumsq == 0 {
    34  		scale = 1
    35  	}
    36  	if scale == 0 {
    37  		scale = 1
    38  		sumsq = 0
    39  	}
    40  
    41  	if n == 0 {
    42  		return scale, sumsq
    43  	}
    44  
    45  	// Compute the sum of squares in 3 accumulators:
    46  	//  - abig: sum of squares scaled down to avoid overflow
    47  	//  - asml: sum of squares scaled up to avoid underflow
    48  	//  - amed: sum of squares that do not require scaling
    49  	// The thresholds and multipliers are:
    50  	//  - values bigger than dtbig are scaled down by dsbig
    51  	//  - values smaller than dtsml are scaled up by dssml
    52  	var (
    53  		isBig            bool
    54  		asml, amed, abig float64
    55  	)
    56  	for i, ix := 0, 0; i < n; i++ {
    57  		ax := math.Abs(x[ix])
    58  		switch {
    59  		case ax > dtbig:
    60  			ax *= dsbig
    61  			abig += ax * ax
    62  			isBig = true
    63  		case ax < dtsml:
    64  			if !isBig {
    65  				ax *= dssml
    66  				asml += ax * ax
    67  			}
    68  		default:
    69  			amed += ax * ax
    70  		}
    71  		ix += incx
    72  	}
    73  	// Put the existing sum of squares into one of the accumulators.
    74  	if sumsq > 0 {
    75  		ax := scale * math.Sqrt(sumsq)
    76  		switch {
    77  		case ax > dtbig:
    78  			// We assume scale >= sqrt( TINY*EPS ) / dsbig, that is, if the
    79  			// scaled sum is big then its scaling factor should not be too
    80  			// small.
    81  			v := scale * dsbig
    82  			abig += (v * v) * sumsq
    83  		case ax < dtsml:
    84  			if !isBig {
    85  				// We assume scale <= sqrt( HUGE ) / dssml, that is, if the
    86  				// scaled sum is small then its scaling factor should not be too
    87  				// big.
    88  				v := scale * dssml
    89  				asml += (v * v) * sumsq
    90  			}
    91  		default:
    92  			amed += ax * ax
    93  		}
    94  	}
    95  	// Combine abig and amed or amed and asml if more than one accumulator was
    96  	// used.
    97  	switch {
    98  	case abig > 0:
    99  		// Combine abig and amed:
   100  		if amed > 0 || math.IsNaN(amed) {
   101  			abig += (amed * dsbig) * dsbig
   102  		}
   103  		scale = 1 / dsbig
   104  		sumsq = abig
   105  	case asml > 0:
   106  		// Combine amed and asml:
   107  		if amed > 0 || math.IsNaN(amed) {
   108  			amed = math.Sqrt(amed)
   109  			asml = math.Sqrt(asml) / dssml
   110  			ymin, ymax := asml, amed
   111  			if asml > amed {
   112  				ymin, ymax = amed, asml
   113  			}
   114  			scale = 1
   115  			sumsq = ymax * ymax * (1 + (ymin/ymax)*(ymin/ymax))
   116  		} else {
   117  			scale = 1 / dssml
   118  			sumsq = asml
   119  		}
   120  	default:
   121  		scale = 1
   122  		sumsq = amed
   123  	}
   124  	return scale, sumsq
   125  }