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

   1 // Copyright ©2017 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/lapack"
  10 )
  11
  12 // Dgerqf computes an RQ factorization of the m×n matrix A,
  13 //  A = R * Q.
  14 // On exit, if m <= n, the upper triangle of the subarray
  15 // A[0:m, n-m:n] contains the m×m upper triangular matrix R.
  16 // If m >= n, the elements on and above the (m-n)-th subdiagonal
  17 // contain the m×n upper trapezoidal matrix R.
  18 // The remaining elements, with tau, represent the
  19 // orthogonal matrix Q as a product of min(m,n) elementary
  20 // reflectors.
  21 //
  22 // The matrix Q is represented as a product of elementary reflectors
  23 //  Q = H_0 H_1 . . . H_{min(m,n)-1}.
  24 // Each H(i) has the form
  25 //  H_i = I - tau_i * v * v^T
  26 // where v is a vector with v[0:n-k+i-1] stored in A[m-k+i, 0:n-k+i-1],
  27 // v[n-k+i:n] = 0 and v[n-k+i] = 1.
  28 //
  29 // tau must have length min(m,n), work must have length max(1, lwork),
  30 // and lwork must be -1 or at least max(1, m), otherwise Dgerqf will panic.
  31 // On exit, work[0] will contain the optimal length for work.
  32 //
  33 // Dgerqf is an internal routine. It is exported for testing purposes.
  34 func (impl Implementation) Dgerqf(m, n int, a []float64, lda int, tau, work []float64, lwork int) {
  35         checkMatrix(m, n, a, lda)
  36
  37         if len(work) < max(1, lwork) {
  38                 panic(shortWork)
  39         }
  40         if lwork != -1 && lwork < max(1, m) {
  41                 panic(badWork)
  42         }
  43
  44         k := min(m, n)
  45         if len(tau) != k {
  46                 panic(badTau)
  47         }
  48
  49         var nb, lwkopt int
  50         if k == 0 {
  51                 lwkopt = 1
  52         } else {
  53                 nb = impl.Ilaenv(1, "DGERQF", " ", m, n, -1, -1)
  54                 lwkopt = m * nb
  55         }
  56         work[0] = float64(lwkopt)
  57
  58         if lwork == -1 {
  59                 return
  60         }
  61
  62         // Return quickly if possible.
  63         if k == 0 {
  64                 return
  65         }
  66
  67         nbmin := 2
  68         nx := 1
  69         iws := m
  70         var ldwork int
  71         if 1 < nb && nb < k {
  72                 // Determine when to cross over from blocked to unblocked code.
  73                 nx = max(0, impl.Ilaenv(3, "DGERQF", " ", m, n, -1, -1))
  74                 if nx < k {
  75                         // Determine whether workspace is large enough for blocked code.
  76                         iws = m * nb
  77                         if lwork < iws {
  78                                 // Not enough workspace to use optimal nb. Reduce
  79                                 // nb and determine the minimum value of nb.
  80                                 nb = lwork / m
  81                                 nbmin = max(2, impl.Ilaenv(2, "DGERQF", " ", m, n, -1, -1))
  82                         }
  83                         ldwork = nb
  84                 }
  85         }
  86
  87         var mu, nu int
  88         if nbmin <= nb && nb < k && nx < k {
  89                 // Use blocked code initially.
  90                 // The last kk rows are handled by the block method.
  91                 ki := ((k - nx - 1) / nb) * nb
  92                 kk := min(k, ki+nb)
  93
  94                 var i int
  95                 for i = k - kk + ki; i >= k-kk; i -= nb {
  96                         ib := min(k-i, nb)
  97
  98                         // Compute the RQ factorization of the current block
  99                         // A[m-k+i:m-k+i+ib-1, 0:n-k+i+ib-1].
 100                         impl.Dgerq2(ib, n-k+i+ib, a[(m-k+i)*lda:], lda, tau[i:], work)
 101                         if m-k+i > 0 {
 102                                 // Form the triangular factor of the block reflector
 103                                 // H = H_{i+ib-1} . . . H_{i+1} H_i.
 104                                 impl.Dlarft(lapack.Backward, lapack.RowWise,
 105                                         n-k+i+ib, ib, a[(m-k+i)*lda:], lda, tau[i:],
 106                                         work, ldwork)
 107
 108                                 // Apply H to A[0:m-k+i-1, 0:n-k+i+ib-1] from the right.
 109                                 impl.Dlarfb(blas.Right, blas.NoTrans, lapack.Backward, lapack.RowWise,
 110                                         m-k+i, n-k+i+ib, ib, a[(m-k+i)*lda:], lda,
 111                                         work, ldwork,
 112                                         a, lda,
 113                                         work[ib*ldwork:], ldwork)
 114                         }
 115                 }
 116                 mu = m - k + i + nb
 117                 nu = n - k + i + nb
 118         } else {
 119                 mu = m
 120                 nu = n
 121         }
 122
 123         // Use unblocked code to factor the last or only block.
 124         if mu > 0 && nu > 0 {
 125                 impl.Dgerq2(mu, nu, a, lda, tau, work)
 126         }
 127         work[0] = float64(iws)
 128 }