vendor/gonum.org/v1/gonum/lapack/gonum/dgebrd.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 // Dgebrd reduces a general m×n matrix A to upper or lower bidiagonal form B by
  13 // an orthogonal transformation:
  14 //  Q^T * A * P = B.
  15 // The diagonal elements of B are stored in d and the off-diagonal elements are stored
  16 // in e. These are additionally stored along the diagonal of A and the off-diagonal
  17 // of A. If m >= n B is an upper-bidiagonal matrix, and if m < n B is a
  18 // lower-bidiagonal matrix.
  19 //
  20 // The remaining elements of A store the data needed to construct Q and P.
  21 // The matrices Q and P are products of elementary reflectors
  22 //  if m >= n, Q = H_0 * H_1 * ... * H_{n-1},
  23 //             P = G_0 * G_1 * ... * G_{n-2},
  24 //  if m < n,  Q = H_0 * H_1 * ... * H_{m-2},
  25 //             P = G_0 * G_1 * ... * G_{m-1},
  26 // where
  27 //  H_i = I - tauQ[i] * v_i * v_i^T,
  28 //  G_i = I - tauP[i] * u_i * u_i^T.
  29 //
  30 // As an example, on exit the entries of A when m = 6, and n = 5
  31 //  [ d   e  u1  u1  u1]
  32 //  [v1   d   e  u2  u2]
  33 //  [v1  v2   d   e  u3]
  34 //  [v1  v2  v3   d   e]
  35 //  [v1  v2  v3  v4   d]
  36 //  [v1  v2  v3  v4  v5]
  37 // and when m = 5, n = 6
  38 //  [ d  u1  u1  u1  u1  u1]
  39 //  [ e   d  u2  u2  u2  u2]
  40 //  [v1   e   d  u3  u3  u3]
  41 //  [v1  v2   e   d  u4  u4]
  42 //  [v1  v2  v3   e   d  u5]
  43 //
  44 // d, tauQ, and tauP must all have length at least min(m,n), and e must have
  45 // length min(m,n) - 1, unless lwork is -1 when there is no check except for
  46 // work which must have a length of at least one.
  47 //
  48 // work is temporary storage, and lwork specifies the usable memory length.
  49 // At minimum, lwork >= max(1,m,n) or be -1 and this function will panic otherwise.
  50 // Dgebrd is blocked decomposition, but the block size is limited
  51 // by the temporary space available. If lwork == -1, instead of performing Dgebrd,
  52 // the optimal work length will be stored into work[0].
  53 //
  54 // Dgebrd is an internal routine. It is exported for testing purposes.
  55 func (impl Implementation) Dgebrd(m, n int, a []float64, lda int, d, e, tauQ, tauP, work []float64, lwork int) {
  56         checkMatrix(m, n, a, lda)
  57         // Calculate optimal work.
  58         nb := impl.Ilaenv(1, "DGEBRD", " ", m, n, -1, -1)
  59         var lworkOpt int
  60         if lwork == -1 {
  61                 if len(work) < 1 {
  62                         panic(badWork)
  63                 }
  64                 lworkOpt = ((m + n) * nb)
  65                 work[0] = float64(max(1, lworkOpt))
  66                 return
  67         }
  68         minmn := min(m, n)
  69         if len(d) < minmn {
  70                 panic(badD)
  71         }
  72         if len(e) < minmn-1 {
  73                 panic(badE)
  74         }
  75         if len(tauQ) < minmn {
  76                 panic(badTauQ)
  77         }
  78         if len(tauP) < minmn {
  79                 panic(badTauP)
  80         }
  81         ws := max(m, n)
  82         if lwork < max(1, ws) {
  83                 panic(badWork)
  84         }
  85         if len(work) < lwork {
  86                 panic(badWork)
  87         }
  88         var nx int
  89         if nb > 1 && nb < minmn {
  90                 nx = max(nb, impl.Ilaenv(3, "DGEBRD", " ", m, n, -1, -1))
  91                 if nx < minmn {
  92                         ws = (m + n) * nb
  93                         if lwork < ws {
  94                                 nbmin := impl.Ilaenv(2, "DGEBRD", " ", m, n, -1, -1)
  95                                 if lwork >= (m+n)*nbmin {
  96                                         nb = lwork / (m + n)
  97                                 } else {
  98                                         nb = minmn
  99                                         nx = minmn
 100                                 }
 101                         }
 102                 }
 103         } else {
 104                 nx = minmn
 105         }
 106         bi := blas64.Implementation()
 107         ldworkx := nb
 108         ldworky := nb
 109         var i int
 110         // Netlib lapack has minmn - nx, but this makes the last nx rows (which by
 111         // default is large) be unblocked. As written here, the blocking is more
 112         // consistent.
 113         for i = 0; i < minmn-nb; i += nb {
 114                 // Reduce rows and columns i:i+nb to bidiagonal form and return
 115                 // the matrices X and Y which are needed to update the unreduced
 116                 // part of the matrix.
 117                 // X is stored in the first m rows of work, y in the next rows.
 118                 x := work[:m*ldworkx]
 119                 y := work[m*ldworkx:]
 120                 impl.Dlabrd(m-i, n-i, nb, a[i*lda+i:], lda,
 121                         d[i:], e[i:], tauQ[i:], tauP[i:],
 122                         x, ldworkx, y, ldworky)
 123
 124                 // Update the trailing submatrix A[i+nb:m,i+nb:n], using an update
 125                 // of the form  A := A - V*Y**T - X*U**T
 126                 bi.Dgemm(blas.NoTrans, blas.Trans, m-i-nb, n-i-nb, nb,
 127                         -1, a[(i+nb)*lda+i:], lda, y[nb*ldworky:], ldworky,
 128                         1, a[(i+nb)*lda+i+nb:], lda)
 129
 130                 bi.Dgemm(blas.NoTrans, blas.NoTrans, m-i-nb, n-i-nb, nb,
 131                         -1, x[nb*ldworkx:], ldworkx, a[i*lda+i+nb:], lda,
 132                         1, a[(i+nb)*lda+i+nb:], lda)
 133
 134                 // Copy diagonal and off-diagonal elements of B back into A.
 135                 if m >= n {
 136                         for j := i; j < i+nb; j++ {
 137                                 a[j*lda+j] = d[j]
 138                                 a[j*lda+j+1] = e[j]
 139                         }
 140                 } else {
 141                         for j := i; j < i+nb; j++ {
 142                                 a[j*lda+j] = d[j]
 143                                 a[(j+1)*lda+j] = e[j]
 144                         }
 145                 }
 146         }
 147         // Use unblocked code to reduce the remainder of the matrix.
 148         impl.Dgebd2(m-i, n-i, a[i*lda+i:], lda, d[i:], e[i:], tauQ[i:], tauP[i:], work)
 149         work[0] = float64(lworkOpt)
 150 }