blob: 3f614b044673bac5f30f26287e7401b6589bb17a [file] [log] [blame]
 // Copyright ©2015 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package gonum import ( "math" "gonum.org/v1/gonum/blas/blas64" ) // Dlarfg generates an elementary reflector for a Householder matrix. It creates // a real elementary reflector of order n such that // H * (alpha) = (beta) // ( x) ( 0) // Hᵀ * H = I // H is represented in the form // H = 1 - tau * (1; v) * (1 vᵀ) // where tau is a real scalar. // // On entry, x contains the vector x, on exit it contains v. // // Dlarfg is an internal routine. It is exported for testing purposes. func (impl Implementation) Dlarfg(n int, alpha float64, x []float64, incX int) (beta, tau float64) { switch { case n < 0: panic(nLT0) case incX <= 0: panic(badIncX) } if n <= 1 { return alpha, 0 } if len(x) < 1+(n-2)*abs(incX) { panic(shortX) } bi := blas64.Implementation() xnorm := bi.Dnrm2(n-1, x, incX) if xnorm == 0 { return alpha, 0 } beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha) safmin := dlamchS / dlamchE knt := 0 if math.Abs(beta) < safmin { // xnorm and beta may be inaccurate, scale x and recompute. rsafmn := 1 / safmin for { knt++ bi.Dscal(n-1, rsafmn, x, incX) beta *= rsafmn alpha *= rsafmn if math.Abs(beta) >= safmin { break } } xnorm = bi.Dnrm2(n-1, x, incX) beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha) } tau = (beta - alpha) / beta bi.Dscal(n-1, 1/(alpha-beta), x, incX) for j := 0; j < knt; j++ { beta *= safmin } return beta, tau }