vendor/gonum.org/v1/gonum/lapack/gonum/dlabrd.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 // Dlabrd reduces the first NB rows and columns of a real general m×n matrix
  13 // A to upper or lower bidiagonal form by an orthogonal transformation
  14 //  Q**T * A * P
  15 // If m >= n, A is reduced to upper bidiagonal form and upon exit the elements
  16 // on and below the diagonal in the first nb columns represent the elementary
  17 // reflectors, and the elements above the diagonal in the first nb rows represent
  18 // the matrix P. If m < n, A is reduced to lower bidiagonal form and the elements
  19 // P is instead stored above the diagonal.
  20 //
  21 // The reduction to bidiagonal form is stored in d and e, where d are the diagonal
  22 // elements, and e are the off-diagonal elements.
  23 //
  24 // The matrices Q and P are products of elementary reflectors
  25 //  Q = H_0 * H_1 * ... * H_{nb-1}
  26 //  P = G_0 * G_1 * ... * G_{nb-1}
  27 // where
  28 //  H_i = I - tauQ[i] * v_i * v_i^T
  29 //  G_i = I - tauP[i] * u_i * u_i^T
  30 //
  31 // As an example, on exit the entries of A when m = 6, n = 5, and nb = 2
  32 //  [ 1   1  u1  u1  u1]
  33 //  [v1   1   1  u2  u2]
  34 //  [v1  v2   a   a   a]
  35 //  [v1  v2   a   a   a]
  36 //  [v1  v2   a   a   a]
  37 //  [v1  v2   a   a   a]
  38 // and when m = 5, n = 6, and nb = 2
  39 //  [ 1  u1  u1  u1  u1  u1]
  40 //  [ 1   1  u2  u2  u2  u2]
  41 //  [v1   1   a   a   a   a]
  42 //  [v1  v2   a   a   a   a]
  43 //  [v1  v2   a   a   a   a]
  44 //
  45 // Dlabrd also returns the matrices X and Y which are used with U and V to
  46 // apply the transformation to the unreduced part of the matrix
  47 //  A := A - V*Y^T - X*U^T
  48 // and returns the matrices X and Y which are needed to apply the
  49 // transformation to the unreduced part of A.
  50 //
  51 // X is an m×nb matrix, Y is an n×nb matrix. d, e, taup, and tauq must all have
  52 // length at least nb. Dlabrd will panic if these size constraints are violated.
  53 //
  54 // Dlabrd is an internal routine. It is exported for testing purposes.
  55 func (impl Implementation) Dlabrd(m, n, nb int, a []float64, lda int, d, e, tauQ, tauP, x []float64, ldx int, y []float64, ldy int) {
  56         checkMatrix(m, n, a, lda)
  57         checkMatrix(m, nb, x, ldx)
  58         checkMatrix(n, nb, y, ldy)
  59         if len(d) < nb {
  60                 panic(badD)
  61         }
  62         if len(e) < nb {
  63                 panic(badE)
  64         }
  65         if len(tauQ) < nb {
  66                 panic(badTauQ)
  67         }
  68         if len(tauP) < nb {
  69                 panic(badTauP)
  70         }
  71         if m <= 0 || n <= 0 {
  72                 return
  73         }
  74         bi := blas64.Implementation()
  75         if m >= n {
  76                 // Reduce to upper bidiagonal form.
  77                 for i := 0; i < nb; i++ {
  78                         bi.Dgemv(blas.NoTrans, m-i, i, -1, a[i*lda:], lda, y[i*ldy:], 1, 1, a[i*lda+i:], lda)
  79                         bi.Dgemv(blas.NoTrans, m-i, i, -1, x[i*ldx:], ldx, a[i:], lda, 1, a[i*lda+i:], lda)
  80
  81                         a[i*lda+i], tauQ[i] = impl.Dlarfg(m-i, a[i*lda+i], a[min(i+1, m-1)*lda+i:], lda)
  82                         d[i] = a[i*lda+i]
  83                         if i < n-1 {
  84                                 // Compute Y[i+1:n, i].
  85                                 a[i*lda+i] = 1
  86                                 bi.Dgemv(blas.Trans, m-i, n-i-1, 1, a[i*lda+i+1:], lda, a[i*lda+i:], lda, 0, y[(i+1)*ldy+i:], ldy)
  87                                 bi.Dgemv(blas.Trans, m-i, i, 1, a[i*lda:], lda, a[i*lda+i:], lda, 0, y[i:], ldy)
  88                                 bi.Dgemv(blas.NoTrans, n-i-1, i, -1, y[(i+1)*ldy:], ldy, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
  89                                 bi.Dgemv(blas.Trans, m-i, i, 1, x[i*ldx:], ldx, a[i*lda+i:], lda, 0, y[i:], ldy)
  90                                 bi.Dgemv(blas.Trans, i, n-i-1, -1, a[i+1:], lda, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
  91                                 bi.Dscal(n-i-1, tauQ[i], y[(i+1)*ldy+i:], ldy)
  92
  93                                 // Update A[i, i+1:n].
  94                                 bi.Dgemv(blas.NoTrans, n-i-1, i+1, -1, y[(i+1)*ldy:], ldy, a[i*lda:], 1, 1, a[i*lda+i+1:], 1)
  95                                 bi.Dgemv(blas.Trans, i, n-i-1, -1, a[i+1:], lda, x[i*ldx:], 1, 1, a[i*lda+i+1:], 1)
  96
  97                                 // Generate reflection P[i] to annihilate A[i, i+2:n].
  98                                 a[i*lda+i+1], tauP[i] = impl.Dlarfg(n-i-1, a[i*lda+i+1], a[i*lda+min(i+2, n-1):], 1)
  99                                 e[i] = a[i*lda+i+1]
 100                                 a[i*lda+i+1] = 1
 101
 102                                 // Compute X[i+1:m, i].
 103                                 bi.Dgemv(blas.NoTrans, m-i-1, n-i-1, 1, a[(i+1)*lda+i+1:], lda, a[i*lda+i+1:], 1, 0, x[(i+1)*ldx+i:], ldx)
 104                                 bi.Dgemv(blas.Trans, n-i-1, i+1, 1, y[(i+1)*ldy:], ldy, a[i*lda+i+1:], 1, 0, x[i:], ldx)
 105                                 bi.Dgemv(blas.NoTrans, m-i-1, i+1, -1, a[(i+1)*lda:], lda, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
 106                                 bi.Dgemv(blas.NoTrans, i, n-i-1, 1, a[i+1:], lda, a[i*lda+i+1:], 1, 0, x[i:], ldx)
 107                                 bi.Dgemv(blas.NoTrans, m-i-1, i, -1, x[(i+1)*ldx:], ldx, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
 108                                 bi.Dscal(m-i-1, tauP[i], x[(i+1)*ldx+i:], ldx)
 109                         }
 110                 }
 111                 return
 112         }
 113         // Reduce to lower bidiagonal form.
 114         for i := 0; i < nb; i++ {
 115                 // Update A[i,i:n]
 116                 bi.Dgemv(blas.NoTrans, n-i, i, -1, y[i*ldy:], ldy, a[i*lda:], 1, 1, a[i*lda+i:], 1)
 117                 bi.Dgemv(blas.Trans, i, n-i, -1, a[i:], lda, x[i*ldx:], 1, 1, a[i*lda+i:], 1)
 118
 119                 // Generate reflection P[i] to annihilate A[i, i+1:n]
 120                 a[i*lda+i], tauP[i] = impl.Dlarfg(n-i, a[i*lda+i], a[i*lda+min(i+1, n-1):], 1)
 121                 d[i] = a[i*lda+i]
 122                 if i < m-1 {
 123                         a[i*lda+i] = 1
 124                         // Compute X[i+1:m, i].
 125                         bi.Dgemv(blas.NoTrans, m-i-1, n-i, 1, a[(i+1)*lda+i:], lda, a[i*lda+i:], 1, 0, x[(i+1)*ldx+i:], ldx)
 126                         bi.Dgemv(blas.Trans, n-i, i, 1, y[i*ldy:], ldy, a[i*lda+i:], 1, 0, x[i:], ldx)
 127                         bi.Dgemv(blas.NoTrans, m-i-1, i, -1, a[(i+1)*lda:], lda, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
 128                         bi.Dgemv(blas.NoTrans, i, n-i, 1, a[i:], lda, a[i*lda+i:], 1, 0, x[i:], ldx)
 129                         bi.Dgemv(blas.NoTrans, m-i-1, i, -1, x[(i+1)*ldx:], ldx, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
 130                         bi.Dscal(m-i-1, tauP[i], x[(i+1)*ldx+i:], ldx)
 131
 132                         // Update A[i+1:m, i].
 133                         bi.Dgemv(blas.NoTrans, m-i-1, i, -1, a[(i+1)*lda:], lda, y[i*ldy:], 1, 1, a[(i+1)*lda+i:], lda)
 134                         bi.Dgemv(blas.NoTrans, m-i-1, i+1, -1, x[(i+1)*ldx:], ldx, a[i:], lda, 1, a[(i+1)*lda+i:], lda)
 135
 136                         // Generate reflection Q[i] to annihilate A[i+2:m, i].
 137                         a[(i+1)*lda+i], tauQ[i] = impl.Dlarfg(m-i-1, a[(i+1)*lda+i], a[min(i+2, m-1)*lda+i:], lda)
 138                         e[i] = a[(i+1)*lda+i]
 139                         a[(i+1)*lda+i] = 1
 140
 141                         // Compute Y[i+1:n, i].
 142                         bi.Dgemv(blas.Trans, m-i-1, n-i-1, 1, a[(i+1)*lda+i+1:], lda, a[(i+1)*lda+i:], lda, 0, y[(i+1)*ldy+i:], ldy)
 143                         bi.Dgemv(blas.Trans, m-i-1, i, 1, a[(i+1)*lda:], lda, a[(i+1)*lda+i:], lda, 0, y[i:], ldy)
 144                         bi.Dgemv(blas.NoTrans, n-i-1, i, -1, y[(i+1)*ldy:], ldy, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
 145                         bi.Dgemv(blas.Trans, m-i-1, i+1, 1, x[(i+1)*ldx:], ldx, a[(i+1)*lda+i:], lda, 0, y[i:], ldy)
 146                         bi.Dgemv(blas.Trans, i+1, n-i-1, -1, a[i+1:], lda, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
 147                         bi.Dscal(n-i-1, tauQ[i], y[(i+1)*ldy+i:], ldy)
 148                 }
 149         }
 150 }