vendor/gonum.org/v1/gonum/lapack/gonum/dormqr.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 // Dormqr multiplies an m×n matrix C by an orthogonal matrix Q as
  13 //  C = Q * C,    if side == blas.Left  and trans == blas.NoTrans,
  14 //  C = Q^T * C,  if side == blas.Left  and trans == blas.Trans,
  15 //  C = C * Q,    if side == blas.Right and trans == blas.NoTrans,
  16 //  C = C * Q^T,  if side == blas.Right and trans == blas.Trans,
  17 // where Q is defined as the product of k elementary reflectors
  18 //  Q = H_0 * H_1 * ... * H_{k-1}.
  19 //
  20 // If side == blas.Left, A is an m×k matrix and 0 <= k <= m.
  21 // If side == blas.Right, A is an n×k matrix and 0 <= k <= n.
  22 // The ith column of A contains the vector which defines the elementary
  23 // reflector H_i and tau[i] contains its scalar factor. tau must have length k
  24 // and Dormqr will panic otherwise. Dgeqrf returns A and tau in the required
  25 // form.
  26 //
  27 // work must have length at least max(1,lwork), and lwork must be at least n if
  28 // side == blas.Left and at least m if side == blas.Right, otherwise Dormqr will
  29 // panic.
  30 //
  31 // work is temporary storage, and lwork specifies the usable memory length. At
  32 // minimum, lwork >= m if side == blas.Left and lwork >= n if side ==
  33 // blas.Right, and this function will panic otherwise. Larger values of lwork
  34 // will generally give better performance. On return, work[0] will contain the
  35 // optimal value of lwork.
  36 //
  37 // If lwork is -1, instead of performing Dormqr, the optimal workspace size will
  38 // be stored into work[0].
  39 func (impl Implementation) Dormqr(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
  40         var nq, nw int
  41         switch side {
  42         default:
  43                 panic(badSide)
  44         case blas.Left:
  45                 nq = m
  46                 nw = n
  47         case blas.Right:
  48                 nq = n
  49                 nw = m
  50         }
  51         switch {
  52         case trans != blas.NoTrans && trans != blas.Trans:
  53                 panic(badTrans)
  54         case m < 0 || n < 0:
  55                 panic(negDimension)
  56         case k < 0 || nq < k:
  57                 panic("lapack: invalid value of k")
  58         case len(work) < lwork:
  59                 panic(shortWork)
  60         case lwork < max(1, nw) && lwork != -1:
  61                 panic(badWork)
  62         }
  63         if lwork != -1 {
  64                 checkMatrix(nq, k, a, lda)
  65                 checkMatrix(m, n, c, ldc)
  66                 if len(tau) != k {
  67                         panic(badTau)
  68                 }
  69         }
  70
  71         if m == 0 || n == 0 || k == 0 {
  72                 work[0] = 1
  73                 return
  74         }
  75
  76         const (
  77                 nbmax = 64
  78                 ldt   = nbmax
  79                 tsize = nbmax * ldt
  80         )
  81         opts := string(side) + string(trans)
  82         nb := min(nbmax, impl.Ilaenv(1, "DORMQR", opts, m, n, k, -1))
  83         lworkopt := max(1, nw)*nb + tsize
  84         if lwork == -1 {
  85                 work[0] = float64(lworkopt)
  86                 return
  87         }
  88
  89         nbmin := 2
  90         if 1 < nb && nb < k {
  91                 if lwork < nw*nb+tsize {
  92                         nb = (lwork - tsize) / nw
  93                         nbmin = max(2, impl.Ilaenv(2, "DORMQR", opts, m, n, k, -1))
  94                 }
  95         }
  96
  97         if nb < nbmin || k <= nb {
  98                 // Call unblocked code.
  99                 impl.Dorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work)
 100                 work[0] = float64(lworkopt)
 101                 return
 102         }
 103
 104         var (
 105                 ldwork = nb
 106                 left   = side == blas.Left
 107                 notran = trans == blas.NoTrans
 108         )
 109         switch {
 110         case left && notran:
 111                 for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
 112                         ib := min(nb, k-i)
 113                         impl.Dlarft(lapack.Forward, lapack.ColumnWise, m-i, ib,
 114                                 a[i*lda+i:], lda,
 115                                 tau[i:],
 116                                 work[:tsize], ldt)
 117                         impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m-i, n, ib,
 118                                 a[i*lda+i:], lda,
 119                                 work[:tsize], ldt,
 120                                 c[i*ldc:], ldc,
 121                                 work[tsize:], ldwork)
 122                 }
 123
 124         case left && !notran:
 125                 for i := 0; i < k; i += nb {
 126                         ib := min(nb, k-i)
 127                         impl.Dlarft(lapack.Forward, lapack.ColumnWise, m-i, ib,
 128                                 a[i*lda+i:], lda,
 129                                 tau[i:],
 130                                 work[:tsize], ldt)
 131                         impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m-i, n, ib,
 132                                 a[i*lda+i:], lda,
 133                                 work[:tsize], ldt,
 134                                 c[i*ldc:], ldc,
 135                                 work[tsize:], ldwork)
 136                 }
 137
 138         case !left && notran:
 139                 for i := 0; i < k; i += nb {
 140                         ib := min(nb, k-i)
 141                         impl.Dlarft(lapack.Forward, lapack.ColumnWise, n-i, ib,
 142                                 a[i*lda+i:], lda,
 143                                 tau[i:],
 144                                 work[:tsize], ldt)
 145                         impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m, n-i, ib,
 146                                 a[i*lda+i:], lda,
 147                                 work[:tsize], ldt,
 148                                 c[i:], ldc,
 149                                 work[tsize:], ldwork)
 150                 }
 151
 152         case !left && !notran:
 153                 for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
 154                         ib := min(nb, k-i)
 155                         impl.Dlarft(lapack.Forward, lapack.ColumnWise, n-i, ib,
 156                                 a[i*lda+i:], lda,
 157                                 tau[i:],
 158                                 work[:tsize], ldt)
 159                         impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m, n-i, ib,
 160                                 a[i*lda+i:], lda,
 161                                 work[:tsize], ldt,
 162                                 c[i:], ldc,
 163                                 work[tsize:], ldwork)
 164                 }
 165         }
 166         work[0] = float64(lworkopt)
 167 }