github.com/gopherd/gonum@v0.0.4/lapack/gonum/lapack.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 "github.com/gopherd/gonum/lapack"
     8  
     9  // Implementation is the native Go implementation of LAPACK routines. It
    10  // is built on top of calls to the return of blas64.Implementation(), so while
    11  // this code is in pure Go, the underlying BLAS implementation may not be.
    12  type Implementation struct{}
    13  
    14  var _ lapack.Float64 = Implementation{}
    15  
    16  func min(a, b int) int {
    17  	if a < b {
    18  		return a
    19  	}
    20  	return b
    21  }
    22  
    23  func max(a, b int) int {
    24  	if a > b {
    25  		return a
    26  	}
    27  	return b
    28  }
    29  
    30  func abs(a int) int {
    31  	if a < 0 {
    32  		return -a
    33  	}
    34  	return a
    35  }
    36  
    37  const (
    38  	// dlamchE is the machine epsilon. For IEEE this is 2^{-53}.
    39  	dlamchE = 0x1p-53
    40  
    41  	// dlamchB is the radix of the machine (the base of the number system).
    42  	dlamchB = 2
    43  
    44  	// dlamchP is base * eps.
    45  	dlamchP = dlamchB * dlamchE
    46  
    47  	// dlamchS is the "safe minimum", that is, the lowest number such that
    48  	// 1/dlamchS does not overflow, or also the smallest normal number.
    49  	// For IEEE this is 2^{-1022}.
    50  	dlamchS = 0x1p-1022
    51  
    52  	// (rtmin,rtmax) is a range of well-scaled numbers whose square
    53  	// or sum of squares is also safe.
    54  	// drtmin is sqrt(dlamchS/dlamchP)
    55  	drtmin = 0x1p-485
    56  	drtmax = 1 / drtmin
    57  
    58  	// Blue's scaling constants
    59  	//
    60  	// An n-vector x is well-scaled if
    61  	//  dtsml ≤ |xᵢ| ≤ dtbig for 0 ≤ i < n and n ≤ 1/dlamchP,
    62  	// where
    63  	//  dtsml = 2^ceil((expmin-1)/2) = 2^ceil((-1021-1)/2) = 2^{-511} = 1.4916681462400413e-154
    64  	//  dtbig = 2^floor((expmax-digits+1)/2) = 2^floor((1024-53+1)/2) = 2^{486} = 1.997919072202235e+146
    65  	// If any xᵢ is not well-scaled, then multiplying small values by dssml and
    66  	// large values by dsbig avoids underflow or overflow when computing the sum
    67  	// of squares \sum_0^{n-1} (xᵢ)².
    68  	//  dssml = 2^{-floor((expmin-digits)/2)} = 2^{-floor((-1021-53)/2)} = 2^537 = 4.4989137945431964e+161
    69  	//  dsbig = 2^{-ceil((expmax+digits-1)/2)} = 2^{-ceil((1024+53-1)/2)} = 2^{-538} = 1.1113793747425387e-162
    70  	//
    71  	// References:
    72  	//  - Anderson E. (2017)
    73  	//    Algorithm 978: Safe Scaling in the Level 1 BLAS
    74  	//    ACM Trans Math Softw 44:1--28
    75  	//    https://doi.org/10.1145/3061665
    76  	//  - Blue, James L. (1978)
    77  	//    A Portable Fortran Program to Find the Euclidean Norm of a Vector
    78  	//    ACM Trans Math Softw 4:15--23
    79  	//    https://doi.org/10.1145/355769.355771
    80  	dtsml = 0x1p-511
    81  	dtbig = 0x1p486
    82  	dssml = 0x1p537
    83  	dsbig = 0x1p-538
    84  )