gonum.org/v1/gonum@v0.14.0/lapack/gonum/dlarfg.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 (
     8  	"math"
     9  
    10  	"gonum.org/v1/gonum/blas/blas64"
    11  )
    12  
    13  // Dlarfg generates an elementary reflector for a Householder matrix. It creates
    14  // a real elementary reflector of order n such that
    15  //
    16  //	H * (alpha) = (beta)
    17  //	    (    x)   (   0)
    18  //	Hᵀ * H = I
    19  //
    20  // H is represented in the form
    21  //
    22  //	H = 1 - tau * (1; v) * (1 vᵀ)
    23  //
    24  // where tau is a real scalar.
    25  //
    26  // On entry, x contains the vector x, on exit it contains v.
    27  //
    28  // Dlarfg is an internal routine. It is exported for testing purposes.
    29  func (impl Implementation) Dlarfg(n int, alpha float64, x []float64, incX int) (beta, tau float64) {
    30  	switch {
    31  	case n < 0:
    32  		panic(nLT0)
    33  	case incX <= 0:
    34  		panic(badIncX)
    35  	}
    36  
    37  	if n <= 1 {
    38  		return alpha, 0
    39  	}
    40  
    41  	if len(x) < 1+(n-2)*abs(incX) {
    42  		panic(shortX)
    43  	}
    44  
    45  	bi := blas64.Implementation()
    46  
    47  	xnorm := bi.Dnrm2(n-1, x, incX)
    48  	if xnorm == 0 {
    49  		return alpha, 0
    50  	}
    51  	beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
    52  	safmin := dlamchS / dlamchE
    53  	knt := 0
    54  	if math.Abs(beta) < safmin {
    55  		// xnorm and beta may be inaccurate, scale x and recompute.
    56  		rsafmn := 1 / safmin
    57  		for {
    58  			knt++
    59  			bi.Dscal(n-1, rsafmn, x, incX)
    60  			beta *= rsafmn
    61  			alpha *= rsafmn
    62  			if math.Abs(beta) >= safmin {
    63  				break
    64  			}
    65  		}
    66  		xnorm = bi.Dnrm2(n-1, x, incX)
    67  		beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
    68  	}
    69  	tau = (beta - alpha) / beta
    70  	bi.Dscal(n-1, 1/(alpha-beta), x, incX)
    71  	for j := 0; j < knt; j++ {
    72  		beta *= safmin
    73  	}
    74  	return beta, tau
    75  }