vendor/gonum.org/v1/gonum/lapack/gonum/dlarfg.go

   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         "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^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 }