vendor/gonum.org/v1/gonum/lapack/gonum/dorgtr.go

   1 // Copyright ©2016 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 // Dorgtr generates a real orthogonal matrix Q which is defined as the product
  10 // of n-1 elementary reflectors of order n as returned by Dsytrd.
  11 //
  12 // The construction of Q depends on the value of uplo:
  13 //  Q = H_{n-1} * ... * H_1 * H_0  if uplo == blas.Upper
  14 //  Q = H_0 * H_1 * ... * H_{n-1}  if uplo == blas.Lower
  15 // where H_i is constructed from the elementary reflectors as computed by Dsytrd.
  16 // See the documentation for Dsytrd for more information.
  17 //
  18 // tau must have length at least n-1, and Dorgtr will panic otherwise.
  19 //
  20 // work is temporary storage, and lwork specifies the usable memory length. At
  21 // minimum, lwork >= max(1,n-1), and Dorgtr will panic otherwise. The amount of blocking
  22 // is limited by the usable length.
  23 // If lwork == -1, instead of computing Dorgtr the optimal work length is stored
  24 // into work[0].
  25 //
  26 // Dorgtr is an internal routine. It is exported for testing purposes.
  27 func (impl Implementation) Dorgtr(uplo blas.Uplo, n int, a []float64, lda int, tau, work []float64, lwork int) {
  28         checkMatrix(n, n, a, lda)
  29         if len(tau) < n-1 {
  30                 panic(badTau)
  31         }
  32         if len(work) < lwork {
  33                 panic(badWork)
  34         }
  35         if lwork < n-1 && lwork != -1 {
  36                 panic(badWork)
  37         }
  38         upper := uplo == blas.Upper
  39         if !upper && uplo != blas.Lower {
  40                 panic(badUplo)
  41         }
  42
  43         if n == 0 {
  44                 work[0] = 1
  45                 return
  46         }
  47
  48         var nb int
  49         if upper {
  50                 nb = impl.Ilaenv(1, "DORGQL", " ", n-1, n-1, n-1, -1)
  51         } else {
  52                 nb = impl.Ilaenv(1, "DORGQR", " ", n-1, n-1, n-1, -1)
  53         }
  54         lworkopt := max(1, n-1) * nb
  55         if lwork == -1 {
  56                 work[0] = float64(lworkopt)
  57                 return
  58         }
  59
  60         if upper {
  61                 // Q was determined by a call to Dsytrd with uplo == blas.Upper.
  62                 // Shift the vectors which define the elementary reflectors one column
  63                 // to the left, and set the last row and column of Q to those of the unit
  64                 // matrix.
  65                 for j := 0; j < n-1; j++ {
  66                         for i := 0; i < j; i++ {
  67                                 a[i*lda+j] = a[i*lda+j+1]
  68                         }
  69                         a[(n-1)*lda+j] = 0
  70                 }
  71                 for i := 0; i < n-1; i++ {
  72                         a[i*lda+n-1] = 0
  73                 }
  74                 a[(n-1)*lda+n-1] = 1
  75
  76                 // Generate Q[0:n-1, 0:n-1].
  77                 impl.Dorgql(n-1, n-1, n-1, a, lda, tau, work, lwork)
  78         } else {
  79                 // Q was determined by a call to Dsytrd with uplo == blas.Upper.
  80                 // Shift the vectors which define the elementary reflectors one column
  81                 // to the right, and set the first row and column of Q to those of the unit
  82                 // matrix.
  83                 for j := n - 1; j > 0; j-- {
  84                         a[j] = 0
  85                         for i := j + 1; i < n; i++ {
  86                                 a[i*lda+j] = a[i*lda+j-1]
  87                         }
  88                 }
  89                 a[0] = 1
  90                 for i := 1; i < n; i++ {
  91                         a[i*lda] = 0
  92                 }
  93                 if n > 1 {
  94                         // Generate Q[1:n, 1:n].
  95                         impl.Dorgqr(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork)
  96                 }
  97         }
  98         work[0] = float64(lworkopt)
  99 }