vendor/gonum.org/v1/gonum/lapack/gonum/dorghr.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 // Dorghr generates an n×n orthogonal matrix Q which is defined as the product
   8 // of ihi-ilo elementary reflectors:
   9 //  Q = H_{ilo} H_{ilo+1} ... H_{ihi-1}.
  10 //
  11 // a and lda represent an n×n matrix that contains the elementary reflectors, as
  12 // returned by Dgehrd. On return, a is overwritten by the n×n orthogonal matrix
  13 // Q. Q will be equal to the identity matrix except in the submatrix
  14 // Q[ilo+1:ihi+1,ilo+1:ihi+1].
  15 //
  16 // ilo and ihi must have the same values as in the previous call of Dgehrd. It
  17 // must hold that
  18 //  0 <= ilo <= ihi < n,  if n > 0,
  19 //  ilo = 0, ihi = -1,    if n == 0.
  20 //
  21 // tau contains the scalar factors of the elementary reflectors, as returned by
  22 // Dgehrd. tau must have length n-1.
  23 //
  24 // work must have length at least max(1,lwork) and lwork must be at least
  25 // ihi-ilo. For optimum performance lwork must be at least (ihi-ilo)*nb where nb
  26 // is the optimal blocksize. On return, work[0] will contain the optimal value
  27 // of lwork.
  28 //
  29 // If lwork == -1, instead of performing Dorghr, only the optimal value of lwork
  30 // will be stored into work[0].
  31 //
  32 // If any requirement on input sizes is not met, Dorghr will panic.
  33 //
  34 // Dorghr is an internal routine. It is exported for testing purposes.
  35 func (impl Implementation) Dorghr(n, ilo, ihi int, a []float64, lda int, tau, work []float64, lwork int) {
  36         checkMatrix(n, n, a, lda)
  37         nh := ihi - ilo
  38         switch {
  39         case ilo < 0 || max(1, n) <= ilo:
  40                 panic(badIlo)
  41         case ihi < min(ilo, n-1) || n <= ihi:
  42                 panic(badIhi)
  43         case lwork < max(1, nh) && lwork != -1:
  44                 panic(badWork)
  45         case len(work) < max(1, lwork):
  46                 panic(shortWork)
  47         }
  48
  49         lwkopt := max(1, nh) * impl.Ilaenv(1, "DORGQR", " ", nh, nh, nh, -1)
  50         if lwork == -1 {
  51                 work[0] = float64(lwkopt)
  52                 return
  53         }
  54
  55         // Quick return if possible.
  56         if n == 0 {
  57                 work[0] = 1
  58                 return
  59         }
  60
  61         // Shift the vectors which define the elementary reflectors one column
  62         // to the right.
  63         for i := ilo + 2; i < ihi+1; i++ {
  64                 copy(a[i*lda+ilo+1:i*lda+i], a[i*lda+ilo:i*lda+i-1])
  65         }
  66         // Set the first ilo+1 and the last n-ihi-1 rows and columns to those of
  67         // the identity matrix.
  68         for i := 0; i < ilo+1; i++ {
  69                 for j := 0; j < n; j++ {
  70                         a[i*lda+j] = 0
  71                 }
  72                 a[i*lda+i] = 1
  73         }
  74         for i := ilo + 1; i < ihi+1; i++ {
  75                 for j := 0; j <= ilo; j++ {
  76                         a[i*lda+j] = 0
  77                 }
  78                 for j := i; j < n; j++ {
  79                         a[i*lda+j] = 0
  80                 }
  81         }
  82         for i := ihi + 1; i < n; i++ {
  83                 for j := 0; j < n; j++ {
  84                         a[i*lda+j] = 0
  85                 }
  86                 a[i*lda+i] = 1
  87         }
  88         if nh > 0 {
  89                 // Generate Q[ilo+1:ihi+1,ilo+1:ihi+1].
  90                 impl.Dorgqr(nh, nh, nh, a[(ilo+1)*lda+ilo+1:], lda, tau[ilo:ihi], work, lwork)
  91         }
  92         work[0] = float64(lwkopt)
  93 }