vendor/gonum.org/v1/gonum/lapack/gonum/dormbr.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 // Dormbr applies a multiplicative update to the matrix C based on a
  13 // decomposition computed by Dgebrd.
  14 //
  15 // Dormbr overwrites the m×n matrix C with
  16 //  Q * C   if vect == lapack.ApplyQ, side == blas.Left, and trans == blas.NoTrans
  17 //  C * Q   if vect == lapack.ApplyQ, side == blas.Right, and trans == blas.NoTrans
  18 //  Q^T * C if vect == lapack.ApplyQ, side == blas.Left, and trans == blas.Trans
  19 //  C * Q^T if vect == lapack.ApplyQ, side == blas.Right, and trans == blas.Trans
  20 //
  21 //  P * C   if vect == lapack.ApplyP, side == blas.Left, and trans == blas.NoTrans
  22 //  C * P   if vect == lapack.ApplyP, side == blas.Right, and trans == blas.NoTrans
  23 //  P^T * C if vect == lapack.ApplyP, side == blas.Left, and trans == blas.Trans
  24 //  C * P^T if vect == lapack.ApplyP, side == blas.Right, and trans == blas.Trans
  25 // where P and Q are the orthogonal matrices determined by Dgebrd when reducing
  26 // a matrix A to bidiagonal form: A = Q * B * P^T. See Dgebrd for the
  27 // definitions of Q and P.
  28 //
  29 // If vect == lapack.ApplyQ, A is assumed to have been an nq×k matrix, while if
  30 // vect == lapack.ApplyP, A is assumed to have been a k×nq matrix. nq = m if
  31 // side == blas.Left, while nq = n if side == blas.Right.
  32 //
  33 // tau must have length min(nq,k), and Dormbr will panic otherwise. tau contains
  34 // the elementary reflectors to construct Q or P depending on the value of
  35 // vect.
  36 //
  37 // work must have length at least max(1,lwork), and lwork must be either -1 or
  38 // at least max(1,n) if side == blas.Left, and at least max(1,m) if side ==
  39 // blas.Right. For optimum performance lwork should be at least n*nb if side ==
  40 // blas.Left, and at least m*nb if side == blas.Right, where nb is the optimal
  41 // block size. On return, work[0] will contain the optimal value of lwork.
  42 //
  43 // If lwork == -1, the function only calculates the optimal value of lwork and
  44 // returns it in work[0].
  45 //
  46 // Dormbr is an internal routine. It is exported for testing purposes.
  47 func (impl Implementation) Dormbr(vect lapack.DecompUpdate, side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
  48         if side != blas.Left && side != blas.Right {
  49                 panic(badSide)
  50         }
  51         if trans != blas.NoTrans && trans != blas.Trans {
  52                 panic(badTrans)
  53         }
  54         if vect != lapack.ApplyP && vect != lapack.ApplyQ {
  55                 panic(badDecompUpdate)
  56         }
  57         nq := n
  58         nw := m
  59         if side == blas.Left {
  60                 nq = m
  61                 nw = n
  62         }
  63         if vect == lapack.ApplyQ {
  64                 checkMatrix(nq, min(nq, k), a, lda)
  65         } else {
  66                 checkMatrix(min(nq, k), nq, a, lda)
  67         }
  68         if len(tau) < min(nq, k) {
  69                 panic(badTau)
  70         }
  71         checkMatrix(m, n, c, ldc)
  72         if len(work) < lwork {
  73                 panic(shortWork)
  74         }
  75         if lwork < max(1, nw) && lwork != -1 {
  76                 panic(badWork)
  77         }
  78
  79         applyQ := vect == lapack.ApplyQ
  80         left := side == blas.Left
  81         var nb int
  82
  83         // The current implementation does not use opts, but a future change may
  84         // use these options so construct them.
  85         var opts string
  86         if side == blas.Left {
  87                 opts = "L"
  88         } else {
  89                 opts = "R"
  90         }
  91         if trans == blas.Trans {
  92                 opts += "T"
  93         } else {
  94                 opts += "N"
  95         }
  96         if applyQ {
  97                 if left {
  98                         nb = impl.Ilaenv(1, "DORMQR", opts, m-1, n, m-1, -1)
  99                 } else {
 100                         nb = impl.Ilaenv(1, "DORMQR", opts, m, n-1, n-1, -1)
 101                 }
 102         } else {
 103                 if left {
 104                         nb = impl.Ilaenv(1, "DORMLQ", opts, m-1, n, m-1, -1)
 105                 } else {
 106                         nb = impl.Ilaenv(1, "DORMLQ", opts, m, n-1, n-1, -1)
 107                 }
 108         }
 109         lworkopt := max(1, nw) * nb
 110         if lwork == -1 {
 111                 work[0] = float64(lworkopt)
 112         }
 113         if applyQ {
 114                 // Change the operation to get Q depending on the size of the initial
 115                 // matrix to Dgebrd. The size matters due to the storage location of
 116                 // the off-diagonal elements.
 117                 if nq >= k {
 118                         impl.Dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork)
 119                 } else if nq > 1 {
 120                         mi := m
 121                         ni := n - 1
 122                         i1 := 0
 123                         i2 := 1
 124                         if left {
 125                                 mi = m - 1
 126                                 ni = n
 127                                 i1 = 1
 128                                 i2 = 0
 129                         }
 130                         impl.Dormqr(side, trans, mi, ni, nq-1, a[1*lda:], lda, tau[:nq-1], c[i1*ldc+i2:], ldc, work, lwork)
 131                 }
 132                 work[0] = float64(lworkopt)
 133                 return
 134         }
 135         transt := blas.Trans
 136         if trans == blas.Trans {
 137                 transt = blas.NoTrans
 138         }
 139         // Change the operation to get P depending on the size of the initial
 140         // matrix to Dgebrd. The size matters due to the storage location of
 141         // the off-diagonal elements.
 142         if nq > k {
 143                 impl.Dormlq(side, transt, m, n, k, a, lda, tau, c, ldc, work, lwork)
 144         } else if nq > 1 {
 145                 mi := m
 146                 ni := n - 1
 147                 i1 := 0
 148                 i2 := 1
 149                 if left {
 150                         mi = m - 1
 151                         ni = n
 152                         i1 = 1
 153                         i2 = 0
 154                 }
 155                 impl.Dormlq(side, transt, mi, ni, nq-1, a[1:], lda, tau, c[i1*ldc+i2:], ldc, work, lwork)
 156         }
 157         work[0] = float64(lworkopt)
 158 }