vendor/gonum.org/v1/gonum/lapack/gonum/dlarf.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         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/blas/blas64"
  10 )
  11
  12 // Dlarf applies an elementary reflector to a general rectangular matrix c.
  13 // This computes
  14 //  c = h * c if side == Left
  15 //  c = c * h if side == right
  16 // where
  17 //  h = 1 - tau * v * v^T
  18 // and c is an m * n matrix.
  19 //
  20 // work is temporary storage of length at least m if side == Left and at least
  21 // n if side == Right. This function will panic if this length requirement is not met.
  22 //
  23 // Dlarf is an internal routine. It is exported for testing purposes.
  24 func (impl Implementation) Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64) {
  25         applyleft := side == blas.Left
  26         if (applyleft && len(work) < n) || (!applyleft && len(work) < m) {
  27                 panic(badWork)
  28         }
  29         checkMatrix(m, n, c, ldc)
  30
  31         // v has length m if applyleft and n otherwise.
  32         lenV := n
  33         if applyleft {
  34                 lenV = m
  35         }
  36
  37         checkVector(lenV, v, incv)
  38
  39         lastv := 0 // last non-zero element of v
  40         lastc := 0 // last non-zero row/column of c
  41         if tau != 0 {
  42                 var i int
  43                 if applyleft {
  44                         lastv = m - 1
  45                 } else {
  46                         lastv = n - 1
  47                 }
  48                 if incv > 0 {
  49                         i = lastv * incv
  50                 }
  51
  52                 // Look for the last non-zero row in v.
  53                 for lastv >= 0 && v[i] == 0 {
  54                         lastv--
  55                         i -= incv
  56                 }
  57                 if applyleft {
  58                         // Scan for the last non-zero column in C[0:lastv, :]
  59                         lastc = impl.Iladlc(lastv+1, n, c, ldc)
  60                 } else {
  61                         // Scan for the last non-zero row in C[:, 0:lastv]
  62                         lastc = impl.Iladlr(m, lastv+1, c, ldc)
  63                 }
  64         }
  65         if lastv == -1 || lastc == -1 {
  66                 return
  67         }
  68         // Sometimes 1-indexing is nicer ...
  69         bi := blas64.Implementation()
  70         if applyleft {
  71                 // Form H * C
  72                 // w[0:lastc+1] = c[1:lastv+1, 1:lastc+1]^T * v[1:lastv+1,1]
  73                 bi.Dgemv(blas.Trans, lastv+1, lastc+1, 1, c, ldc, v, incv, 0, work, 1)
  74                 // c[0: lastv, 0: lastc] = c[...] - w[0:lastv, 1] * v[1:lastc, 1]^T
  75                 bi.Dger(lastv+1, lastc+1, -tau, v, incv, work, 1, c, ldc)
  76                 return
  77         }
  78         // Form C*H
  79         // w[0:lastc+1,1] := c[0:lastc+1,0:lastv+1] * v[0:lastv+1,1]
  80         bi.Dgemv(blas.NoTrans, lastc+1, lastv+1, 1, c, ldc, v, incv, 0, work, 1)
  81         // c[0:lastc+1,0:lastv+1] = c[...] - w[0:lastc+1,0] * v[0:lastv+1,0]^T
  82         bi.Dger(lastc+1, lastv+1, -tau, work, 1, v, incv, c, ldc)
  83 }