vendor/gonum.org/v1/gonum/lapack/gonum/dorgbr.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/lapack"
   8
   9 // Dorgbr generates one of the matrices Q or P^T computed by Dgebrd
  10 // computed from the decomposition Dgebrd. See Dgebd2 for the description of
  11 // Q and P^T.
  12 //
  13 // If vect == lapack.ApplyQ, then a is assumed to have been an m×k matrix and
  14 // Q is of order m. If m >= k, then Dorgbr returns the first n columns of Q
  15 // where m >= n >= k. If m < k, then Dorgbr returns Q as an m×m matrix.
  16 //
  17 // If vect == lapack.ApplyP, then A is assumed to have been a k×n matrix, and
  18 // P^T is of order n. If k < n, then Dorgbr returns the first m rows of P^T,
  19 // where n >= m >= k. If k >= n, then Dorgbr returns P^T as an n×n matrix.
  20 //
  21 // Dorgbr is an internal routine. It is exported for testing purposes.
  22 func (impl Implementation) Dorgbr(vect lapack.DecompUpdate, m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
  23         mn := min(m, n)
  24         wantq := vect == lapack.ApplyQ
  25         if wantq {
  26                 if m < n || n < min(m, k) || m < min(m, k) {
  27                         panic(badDims)
  28                 }
  29         } else {
  30                 if n < m || m < min(n, k) || n < min(n, k) {
  31                         panic(badDims)
  32                 }
  33         }
  34         if wantq {
  35                 if m >= k {
  36                         checkMatrix(m, k, a, lda)
  37                 } else {
  38                         checkMatrix(m, m, a, lda)
  39                 }
  40         } else {
  41                 if n >= k {
  42                         checkMatrix(k, n, a, lda)
  43                 } else {
  44                         checkMatrix(n, n, a, lda)
  45                 }
  46         }
  47         work[0] = 1
  48         if wantq {
  49                 if m >= k {
  50                         impl.Dorgqr(m, n, k, a, lda, tau, work, -1)
  51                 } else if m > 1 {
  52                         impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, -1)
  53                 }
  54         } else {
  55                 if k < n {
  56                         impl.Dorglq(m, n, k, a, lda, tau, work, -1)
  57                 } else if n > 1 {
  58                         impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, -1)
  59                 }
  60         }
  61         lworkopt := int(work[0])
  62         lworkopt = max(lworkopt, mn)
  63         if lwork == -1 {
  64                 work[0] = float64(lworkopt)
  65                 return
  66         }
  67         if len(work) < lwork {
  68                 panic(badWork)
  69         }
  70         if lwork < mn {
  71                 panic(badWork)
  72         }
  73         if m == 0 || n == 0 {
  74                 work[0] = 1
  75                 return
  76         }
  77         if wantq {
  78                 // Form Q, determined by a call to Dgebrd to reduce an m×k matrix.
  79                 if m >= k {
  80                         impl.Dorgqr(m, n, k, a, lda, tau, work, lwork)
  81                 } else {
  82                         // Shift the vectors which define the elementary reflectors one
  83                         // column to the right, and set the first row and column of Q to
  84                         // those of the unit matrix.
  85                         for j := m - 1; j >= 1; j-- {
  86                                 a[j] = 0
  87                                 for i := j + 1; i < m; i++ {
  88                                         a[i*lda+j] = a[i*lda+j-1]
  89                                 }
  90                         }
  91                         a[0] = 1
  92                         for i := 1; i < m; i++ {
  93                                 a[i*lda] = 0
  94                         }
  95                         if m > 1 {
  96                                 // Form Q[1:m-1, 1:m-1]
  97                                 impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, lwork)
  98                         }
  99                 }
 100         } else {
 101                 // Form P^T, determined by a call to Dgebrd to reduce a k×n matrix.
 102                 if k < n {
 103                         impl.Dorglq(m, n, k, a, lda, tau, work, lwork)
 104                 } else {
 105                         // Shift the vectors which define the elementary reflectors one
 106                         // row downward, and set the first row and column of P^T to
 107                         // those of the unit matrix.
 108                         a[0] = 1
 109                         for i := 1; i < n; i++ {
 110                                 a[i*lda] = 0
 111                         }
 112                         for j := 1; j < n; j++ {
 113                                 for i := j - 1; i >= 1; i-- {
 114                                         a[i*lda+j] = a[(i-1)*lda+j]
 115                                 }
 116                                 a[j] = 0
 117                         }
 118                         if n > 1 {
 119                                 impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork)
 120                         }
 121                 }
 122         }
 123         work[0] = float64(lworkopt)
 124 }