vendor/gonum.org/v1/gonum/lapack/gonum/dorglq.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/lapack"
  10 )
  11
  12 // Dorglq generates an m×n matrix Q with orthonormal columns defined by the
  13 // product of elementary reflectors as computed by Dgelqf.
  14 //  Q = H_0 * H_1 * ... * H_{k-1}
  15 // Dorglq is the blocked version of Dorgl2 that makes greater use of level-3 BLAS
  16 // routines.
  17 //
  18 // len(tau) >= k, 0 <= k <= n, and 0 <= n <= m.
  19 //
  20 // work is temporary storage, and lwork specifies the usable memory length. At minimum,
  21 // lwork >= m, and the amount of blocking is limited by the usable length.
  22 // If lwork == -1, instead of computing Dorglq the optimal work length is stored
  23 // into work[0].
  24 //
  25 // Dorglq will panic if the conditions on input values are not met.
  26 //
  27 // Dorglq is an internal routine. It is exported for testing purposes.
  28 func (impl Implementation) Dorglq(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
  29         nb := impl.Ilaenv(1, "DORGLQ", " ", m, n, k, -1)
  30         // work is treated as an n×nb matrix
  31         if lwork == -1 {
  32                 work[0] = float64(max(1, m) * nb)
  33                 return
  34         }
  35         checkMatrix(m, n, a, lda)
  36         if k < 0 {
  37                 panic(kLT0)
  38         }
  39         if k > m {
  40                 panic(kGTM)
  41         }
  42         if m > n {
  43                 panic(nLTM)
  44         }
  45         if len(tau) < k {
  46                 panic(badTau)
  47         }
  48         if len(work) < lwork {
  49                 panic(shortWork)
  50         }
  51         if lwork < m {
  52                 panic(badWork)
  53         }
  54         if m == 0 {
  55                 return
  56         }
  57         nbmin := 2 // Minimum number of blocks
  58         var nx int // Minimum number of rows
  59         iws := m   // Length of work needed
  60         var ldwork int
  61         if nb > 1 && nb < k {
  62                 nx = max(0, impl.Ilaenv(3, "DORGLQ", " ", m, n, k, -1))
  63                 if nx < k {
  64                         ldwork = nb
  65                         iws = m * ldwork
  66                         if lwork < iws {
  67                                 nb = lwork / m
  68                                 ldwork = nb
  69                                 nbmin = max(2, impl.Ilaenv(2, "DORGLQ", " ", m, n, k, -1))
  70                         }
  71                 }
  72         }
  73         var ki, kk int
  74         if nb >= nbmin && nb < k && nx < k {
  75                 // The first kk rows are handled by the blocked method.
  76                 // Note: lapack has nx here, but this means the last nx rows are handled
  77                 // serially which could be quite different than nb.
  78                 ki = ((k - nb - 1) / nb) * nb
  79                 kk = min(k, ki+nb)
  80                 for i := kk; i < m; i++ {
  81                         for j := 0; j < kk; j++ {
  82                                 a[i*lda+j] = 0
  83                         }
  84                 }
  85         }
  86         if kk < m {
  87                 // Perform the operation on colums kk to the end.
  88                 impl.Dorgl2(m-kk, n-kk, k-kk, a[kk*lda+kk:], lda, tau[kk:], work)
  89         }
  90         if kk == 0 {
  91                 return
  92         }
  93         // Perform the operation on column-blocks
  94         for i := ki; i >= 0; i -= nb {
  95                 ib := min(nb, k-i)
  96                 if i+ib < m {
  97                         impl.Dlarft(lapack.Forward, lapack.RowWise,
  98                                 n-i, ib,
  99                                 a[i*lda+i:], lda,
 100                                 tau[i:],
 101                                 work, ldwork)
 102
 103                         impl.Dlarfb(blas.Right, blas.Trans, lapack.Forward, lapack.RowWise,
 104                                 m-i-ib, n-i, ib,
 105                                 a[i*lda+i:], lda,
 106                                 work, ldwork,
 107                                 a[(i+ib)*lda+i:], lda,
 108                                 work[ib*ldwork:], ldwork)
 109                 }
 110                 impl.Dorgl2(ib, n-i, ib, a[i*lda+i:], lda, tau[i:], work)
 111                 for l := i; l < i+ib; l++ {
 112                         for j := 0; j < i; j++ {
 113                                 a[l*lda+j] = 0
 114                         }
 115                 }
 116         }
 117 }