gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/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  			if scale > 1 {
    81  				scale *= dsbig
    82  				abig += scale * (scale * sumsq)
    83  			} else {
    84  				// sumsq > dtbig^2 => (dsbig * (dsbig * sumsq)) is representable.
    85  				abig += scale * (scale * (dsbig * (dsbig * sumsq)))
    86  			}
    87  		case ax < dtsml:
    88  			if !isBig {
    89  				if scale < 1 {
    90  					scale *= dssml
    91  					asml += scale * (scale * sumsq)
    92  				} else {
    93  					// sumsq < dtsml^2 => (dssml * (dssml * sumsq)) is representable.
    94  					asml += scale * (scale * (dssml * (dssml * sumsq)))
    95  				}
    96  			}
    97  		default:
    98  			amed += scale * (scale * sumsq)
    99  		}
   100  	}
   101  	// Combine abig and amed or amed and asml if more than one accumulator was
   102  	// used.
   103  	switch {
   104  	case abig > 0:
   105  		// Combine abig and amed:
   106  		if amed > 0 || math.IsNaN(amed) {
   107  			abig += (amed * dsbig) * dsbig
   108  		}
   109  		scale = 1 / dsbig
   110  		sumsq = abig
   111  	case asml > 0:
   112  		// Combine amed and asml:
   113  		if amed > 0 || math.IsNaN(amed) {
   114  			amed = math.Sqrt(amed)
   115  			asml = math.Sqrt(asml) / dssml
   116  			ymin, ymax := asml, amed
   117  			if asml > amed {
   118  				ymin, ymax = amed, asml
   119  			}
   120  			scale = 1
   121  			sumsq = ymax * ymax * (1 + (ymin/ymax)*(ymin/ymax))
   122  		} else {
   123  			scale = 1 / dssml
   124  			sumsq = asml
   125  		}
   126  	default:
   127  		scale = 1
   128  		sumsq = amed
   129  	}
   130  	return scale, sumsq
   131  }