github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/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^T * H = I
    18  // H is represented in the form
    19  //  H = 1 - tau * (1; v) * (1 v^T)
    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  	if n < 0 {
    27  		panic(nLT0)
    28  	}
    29  	if n <= 1 {
    30  		return alpha, 0
    31  	}
    32  	checkVector(n-1, x, incX)
    33  	bi := blas64.Implementation()
    34  	xnorm := bi.Dnrm2(n-1, x, incX)
    35  	if xnorm == 0 {
    36  		return alpha, 0
    37  	}
    38  	beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
    39  	safmin := dlamchS / dlamchE
    40  	knt := 0
    41  	if math.Abs(beta) < safmin {
    42  		// xnorm and beta may be inaccurate, scale x and recompute.
    43  		rsafmn := 1 / safmin
    44  		for {
    45  			knt++
    46  			bi.Dscal(n-1, rsafmn, x, incX)
    47  			beta *= rsafmn
    48  			alpha *= rsafmn
    49  			if math.Abs(beta) >= safmin {
    50  				break
    51  			}
    52  		}
    53  		xnorm = bi.Dnrm2(n-1, x, incX)
    54  		beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
    55  	}
    56  	tau = (beta - alpha) / beta
    57  	bi.Dscal(n-1, 1/(alpha-beta), x, incX)
    58  	for j := 0; j < knt; j++ {
    59  		beta *= safmin
    60  	}
    61  	return beta, tau
    62  }