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