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

   1 // Copyright ©2016 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 // Dorgql generates the m×n matrix Q with orthonormal columns defined as the
  13 // last n columns of a product of k elementary reflectors of order m
  14 //  Q = H_{k-1} * ... * H_1 * H_0.
  15 //
  16 // It must hold that
  17 //  0 <= k <= n <= m,
  18 // and Dorgql will panic otherwise.
  19 //
  20 // On entry, the (n-k+i)-th column of A must contain the vector which defines
  21 // the elementary reflector H_i, for i=0,...,k-1, and tau[i] must contain its
  22 // scalar factor. On return, a contains the m×n matrix Q.
  23 //
  24 // tau must have length at least k, and Dorgql will panic otherwise.
  25 //
  26 // work must have length at least max(1,lwork), and lwork must be at least
  27 // max(1,n), otherwise Dorgql will panic. For optimum performance lwork must
  28 // be a sufficiently large multiple of n.
  29 //
  30 // If lwork == -1, instead of computing Dorgql the optimal work length is stored
  31 // into work[0].
  32 //
  33 // Dorgql is an internal routine. It is exported for testing purposes.
  34 func (impl Implementation) Dorgql(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
  35         switch {
  36         case n < 0:
  37                 panic(nLT0)
  38         case m < n:
  39                 panic(mLTN)
  40         case k < 0:
  41                 panic(kLT0)
  42         case k > n:
  43                 panic(kGTN)
  44         case lwork < max(1, n) && lwork != -1:
  45                 panic(badWork)
  46         case len(work) < lwork:
  47                 panic(shortWork)
  48         }
  49         if lwork != -1 {
  50                 checkMatrix(m, n, a, lda)
  51                 if len(tau) < k {
  52                         panic(badTau)
  53                 }
  54         }
  55
  56         if n == 0 {
  57                 work[0] = 1
  58                 return
  59         }
  60
  61         nb := impl.Ilaenv(1, "DORGQL", " ", m, n, k, -1)
  62         if lwork == -1 {
  63                 work[0] = float64(n * nb)
  64                 return
  65         }
  66
  67         nbmin := 2
  68         var nx, ldwork int
  69         iws := n
  70         if nb > 1 && nb < k {
  71                 // Determine when to cross over from blocked to unblocked code.
  72                 nx = max(0, impl.Ilaenv(3, "DORGQL", " ", m, n, k, -1))
  73                 if nx < k {
  74                         // Determine if workspace is large enough for blocked code.
  75                         iws = n * nb
  76                         if lwork < iws {
  77                                 // Not enough workspace to use optimal nb: reduce nb and determine
  78                                 // the minimum value of nb.
  79                                 nb = lwork / n
  80                                 nbmin = max(2, impl.Ilaenv(2, "DORGQL", " ", m, n, k, -1))
  81                         }
  82                         ldwork = nb
  83                 }
  84         }
  85
  86         var kk int
  87         if nb >= nbmin && nb < k && nx < k {
  88                 // Use blocked code after the first block. The last kk columns are handled
  89                 // by the block method.
  90                 kk = min(k, ((k-nx+nb-1)/nb)*nb)
  91
  92                 // Set A(m-kk:m, 0:n-kk) to zero.
  93                 for i := m - kk; i < m; i++ {
  94                         for j := 0; j < n-kk; j++ {
  95                                 a[i*lda+j] = 0
  96                         }
  97                 }
  98         }
  99
 100         // Use unblocked code for the first or only block.
 101         impl.Dorg2l(m-kk, n-kk, k-kk, a, lda, tau, work)
 102         if kk > 0 {
 103                 // Use blocked code.
 104                 for i := k - kk; i < k; i += nb {
 105                         ib := min(nb, k-i)
 106                         if n-k+i > 0 {
 107                                 // Form the triangular factor of the block reflector
 108                                 // H = H_{i+ib-1} * ... * H_{i+1} * H_i.
 109                                 impl.Dlarft(lapack.Backward, lapack.ColumnWise, m-k+i+ib, ib,
 110                                         a[n-k+i:], lda, tau[i:], work, ldwork)
 111
 112                                 // Apply H to A[0:m-k+i+ib, 0:n-k+i] from the left.
 113                                 impl.Dlarfb(blas.Left, blas.NoTrans, lapack.Backward, lapack.ColumnWise,
 114                                         m-k+i+ib, n-k+i, ib, a[n-k+i:], lda, work, ldwork,
 115                                         a, lda, work[ib*ldwork:], ldwork)
 116                         }
 117
 118                         // Apply H to rows 0:m-k+i+ib of current block.
 119                         impl.Dorg2l(m-k+i+ib, ib, ib, a[n-k+i:], lda, tau[i:], work)
 120
 121                         // Set rows m-k+i+ib:m of current block to zero.
 122                         for j := n - k + i; j < n-k+i+ib; j++ {
 123                                 for l := m - k + i + ib; l < m; l++ {
 124                                         a[l*lda+j] = 0
 125                                 }
 126                         }
 127                 }
 128         }
 129         work[0] = float64(iws)
 130 }