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
5package gonum
6
7import (
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//  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.
25func (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}
72