vendor/gonum.org/v1/gonum/lapack/gonum/dormlq.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 // Dormlq multiplies the matrix C by the orthogonal matrix Q defined by the
  13 // slices a and tau. A and tau are as returned from Dgelqf.
  14 //  C = Q * C    if side == blas.Left and trans == blas.NoTrans
  15 //  C = Q^T * C  if side == blas.Left and trans == blas.Trans
  16 //  C = C * Q    if side == blas.Right and trans == blas.NoTrans
  17 //  C = C * Q^T  if side == blas.Right and trans == blas.Trans
  18 // If side == blas.Left, A is a matrix of side k×m, and if side == blas.Right
  19 // A is of size k×n. This uses a blocked algorithm.
  20 //
  21 // work is temporary storage, and lwork specifies the usable memory length.
  22 // At minimum, lwork >= m if side == blas.Left and lwork >= n if side == blas.Right,
  23 // and this function will panic otherwise.
  24 // Dormlq uses a block algorithm, but the block size is limited
  25 // by the temporary space available. If lwork == -1, instead of performing Dormlq,
  26 // the optimal work length will be stored into work[0].
  27 //
  28 // tau contains the Householder scales and must have length at least k, and
  29 // this function will panic otherwise.
  30 func (impl Implementation) Dormlq(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
  31         if side != blas.Left && side != blas.Right {
  32                 panic(badSide)
  33         }
  34         if trans != blas.Trans && trans != blas.NoTrans {
  35                 panic(badTrans)
  36         }
  37         left := side == blas.Left
  38         if left {
  39                 checkMatrix(k, m, a, lda)
  40         } else {
  41                 checkMatrix(k, n, a, lda)
  42         }
  43         checkMatrix(m, n, c, ldc)
  44         if len(tau) < k {
  45                 panic(badTau)
  46         }
  47         if len(work) < lwork {
  48                 panic(shortWork)
  49         }
  50         nw := m
  51         if left {
  52                 nw = n
  53         }
  54         if lwork < max(1, nw) && lwork != -1 {
  55                 panic(badWork)
  56         }
  57
  58         if m == 0 || n == 0 || k == 0 {
  59                 work[0] = 1
  60                 return
  61         }
  62
  63         const (
  64                 nbmax = 64
  65                 ldt   = nbmax
  66                 tsize = nbmax * ldt
  67         )
  68         opts := string(side) + string(trans)
  69         nb := min(nbmax, impl.Ilaenv(1, "DORMLQ", opts, m, n, k, -1))
  70         lworkopt := max(1, nw)*nb + tsize
  71         if lwork == -1 {
  72                 work[0] = float64(lworkopt)
  73                 return
  74         }
  75
  76         nbmin := 2
  77         if 1 < nb && nb < k {
  78                 iws := nw*nb + tsize
  79                 if lwork < iws {
  80                         nb = (lwork - tsize) / nw
  81                         nbmin = max(2, impl.Ilaenv(2, "DORMLQ", opts, m, n, k, -1))
  82                 }
  83         }
  84         if nb < nbmin || k <= nb {
  85                 // Call unblocked code.
  86                 impl.Dorml2(side, trans, m, n, k, a, lda, tau, c, ldc, work)
  87                 work[0] = float64(lworkopt)
  88                 return
  89         }
  90
  91         t := work[:tsize]
  92         wrk := work[tsize:]
  93         ldwrk := nb
  94
  95         notran := trans == blas.NoTrans
  96         transt := blas.NoTrans
  97         if notran {
  98                 transt = blas.Trans
  99         }
 100
 101         switch {
 102         case left && notran:
 103                 for i := 0; i < k; i += nb {
 104                         ib := min(nb, k-i)
 105                         impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
 106                                 a[i*lda+i:], lda,
 107                                 tau[i:],
 108                                 t, ldt)
 109                         impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
 110                                 a[i*lda+i:], lda,
 111                                 t, ldt,
 112                                 c[i*ldc:], ldc,
 113                                 wrk, ldwrk)
 114                 }
 115
 116         case left && !notran:
 117                 for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
 118                         ib := min(nb, k-i)
 119                         impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
 120                                 a[i*lda+i:], lda,
 121                                 tau[i:],
 122                                 t, ldt)
 123                         impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
 124                                 a[i*lda+i:], lda,
 125                                 t, ldt,
 126                                 c[i*ldc:], ldc,
 127                                 wrk, ldwrk)
 128                 }
 129
 130         case !left && notran:
 131                 for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
 132                         ib := min(nb, k-i)
 133                         impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
 134                                 a[i*lda+i:], lda,
 135                                 tau[i:],
 136                                 t, ldt)
 137                         impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
 138                                 a[i*lda+i:], lda,
 139                                 t, ldt,
 140                                 c[i:], ldc,
 141                                 wrk, ldwrk)
 142                 }
 143
 144         case !left && !notran:
 145                 for i := 0; i < k; i += nb {
 146                         ib := min(nb, k-i)
 147                         impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
 148                                 a[i*lda+i:], lda,
 149                                 tau[i:],
 150                                 t, ldt)
 151                         impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
 152                                 a[i*lda+i:], lda,
 153                                 t, ldt,
 154                                 c[i:], ldc,
 155                                 wrk, ldwrk)
 156                 }
 157         }
 158         work[0] = float64(lworkopt)
 159 }