--- /dev/null
+// Copyright ©2015 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package gonum
+
+import "gonum.org/v1/gonum/lapack"
+
+// Dorgbr generates one of the matrices Q or P^T computed by Dgebrd
+// computed from the decomposition Dgebrd. See Dgebd2 for the description of
+// Q and P^T.
+//
+// If vect == lapack.ApplyQ, then a is assumed to have been an m×k matrix and
+// Q is of order m. If m >= k, then Dorgbr returns the first n columns of Q
+// where m >= n >= k. If m < k, then Dorgbr returns Q as an m×m matrix.
+//
+// If vect == lapack.ApplyP, then A is assumed to have been a k×n matrix, and
+// P^T is of order n. If k < n, then Dorgbr returns the first m rows of P^T,
+// where n >= m >= k. If k >= n, then Dorgbr returns P^T as an n×n matrix.
+//
+// Dorgbr is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dorgbr(vect lapack.DecompUpdate, m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
+ mn := min(m, n)
+ wantq := vect == lapack.ApplyQ
+ if wantq {
+ if m < n || n < min(m, k) || m < min(m, k) {
+ panic(badDims)
+ }
+ } else {
+ if n < m || m < min(n, k) || n < min(n, k) {
+ panic(badDims)
+ }
+ }
+ if wantq {
+ if m >= k {
+ checkMatrix(m, k, a, lda)
+ } else {
+ checkMatrix(m, m, a, lda)
+ }
+ } else {
+ if n >= k {
+ checkMatrix(k, n, a, lda)
+ } else {
+ checkMatrix(n, n, a, lda)
+ }
+ }
+ work[0] = 1
+ if wantq {
+ if m >= k {
+ impl.Dorgqr(m, n, k, a, lda, tau, work, -1)
+ } else if m > 1 {
+ impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, -1)
+ }
+ } else {
+ if k < n {
+ impl.Dorglq(m, n, k, a, lda, tau, work, -1)
+ } else if n > 1 {
+ impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, -1)
+ }
+ }
+ lworkopt := int(work[0])
+ lworkopt = max(lworkopt, mn)
+ if lwork == -1 {
+ work[0] = float64(lworkopt)
+ return
+ }
+ if len(work) < lwork {
+ panic(badWork)
+ }
+ if lwork < mn {
+ panic(badWork)
+ }
+ if m == 0 || n == 0 {
+ work[0] = 1
+ return
+ }
+ if wantq {
+ // Form Q, determined by a call to Dgebrd to reduce an m×k matrix.
+ if m >= k {
+ impl.Dorgqr(m, n, k, a, lda, tau, work, lwork)
+ } else {
+ // Shift the vectors which define the elementary reflectors one
+ // column to the right, and set the first row and column of Q to
+ // those of the unit matrix.
+ for j := m - 1; j >= 1; j-- {
+ a[j] = 0
+ for i := j + 1; i < m; i++ {
+ a[i*lda+j] = a[i*lda+j-1]
+ }
+ }
+ a[0] = 1
+ for i := 1; i < m; i++ {
+ a[i*lda] = 0
+ }
+ if m > 1 {
+ // Form Q[1:m-1, 1:m-1]
+ impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, lwork)
+ }
+ }
+ } else {
+ // Form P^T, determined by a call to Dgebrd to reduce a k×n matrix.
+ if k < n {
+ impl.Dorglq(m, n, k, a, lda, tau, work, lwork)
+ } else {
+ // Shift the vectors which define the elementary reflectors one
+ // row downward, and set the first row and column of P^T to
+ // those of the unit matrix.
+ a[0] = 1
+ for i := 1; i < n; i++ {
+ a[i*lda] = 0
+ }
+ for j := 1; j < n; j++ {
+ for i := j - 1; i >= 1; i-- {
+ a[i*lda+j] = a[(i-1)*lda+j]
+ }
+ a[j] = 0
+ }
+ if n > 1 {
+ impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork)
+ }
+ }
+ }
+ work[0] = float64(lworkopt)
+}
OSDN Git Service
