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