1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
// 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
}
|
