vendor/gonum.org/v1/gonum/lapack/gonum/dormr2.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 "gonum.org/v1/gonum/blas"
   8
   9 // Dormr2 multiplies a general matrix C by an orthogonal matrix from a RQ factorization
  10 // determined by Dgerqf.
  11 //  C = Q * C    if side == blas.Left and trans == blas.NoTrans
  12 //  C = Q^T * C  if side == blas.Left and trans == blas.Trans
  13 //  C = C * Q    if side == blas.Right and trans == blas.NoTrans
  14 //  C = C * Q^T  if side == blas.Right and trans == blas.Trans
  15 // If side == blas.Left, a is a matrix of size k×m, and if side == blas.Right
  16 // a is of size k×n.
  17 //
  18 // tau contains the Householder factors and is of length at least k and this function
  19 // will panic otherwise.
  20 //
  21 // work is temporary storage of length at least n if side == blas.Left
  22 // and at least m if side == blas.Right and this function will panic otherwise.
  23 //
  24 // Dormr2 is an internal routine. It is exported for testing purposes.
  25 func (impl Implementation) Dormr2(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64) {
  26         if side != blas.Left && side != blas.Right {
  27                 panic(badSide)
  28         }
  29         if trans != blas.Trans && trans != blas.NoTrans {
  30                 panic(badTrans)
  31         }
  32
  33         left := side == blas.Left
  34         notran := trans == blas.NoTrans
  35         if left {
  36                 if k > m {
  37                         panic(kGTM)
  38                 }
  39                 checkMatrix(k, m, a, lda)
  40                 if len(work) < n {
  41                         panic(badWork)
  42                 }
  43         } else {
  44                 if k > n {
  45                         panic(kGTN)
  46                 }
  47                 checkMatrix(k, n, a, lda)
  48                 if len(work) < m {
  49                         panic(badWork)
  50                 }
  51         }
  52         if len(tau) < k {
  53                 panic(badTau)
  54         }
  55         checkMatrix(m, n, c, ldc)
  56
  57         if m == 0 || n == 0 || k == 0 {
  58                 return
  59         }
  60         if left {
  61                 if notran {
  62                         for i := k - 1; i >= 0; i-- {
  63                                 aii := a[i*lda+(m-k+i)]
  64                                 a[i*lda+(m-k+i)] = 1
  65                                 impl.Dlarf(side, m-k+i+1, n, a[i*lda:], 1, tau[i], c, ldc, work)
  66                                 a[i*lda+(m-k+i)] = aii
  67                         }
  68                         return
  69                 }
  70                 for i := 0; i < k; i++ {
  71                         aii := a[i*lda+(m-k+i)]
  72                         a[i*lda+(m-k+i)] = 1
  73                         impl.Dlarf(side, m-k+i+1, n, a[i*lda:], 1, tau[i], c, ldc, work)
  74                         a[i*lda+(m-k+i)] = aii
  75                 }
  76                 return
  77         }
  78         if notran {
  79                 for i := 0; i < k; i++ {
  80                         aii := a[i*lda+(n-k+i)]
  81                         a[i*lda+(n-k+i)] = 1
  82                         impl.Dlarf(side, m, n-k+i+1, a[i*lda:], 1, tau[i], c, ldc, work)
  83                         a[i*lda+(n-k+i)] = aii
  84                 }
  85                 return
  86         }
  87         for i := k - 1; i >= 0; i-- {
  88                 aii := a[i*lda+(n-k+i)]
  89                 a[i*lda+(n-k+i)] = 1
  90                 impl.Dlarf(side, m, n-k+i+1, a[i*lda:], 1, tau[i], c, ldc, work)
  91                 a[i*lda+(n-k+i)] = aii
  92         }
  93 }