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.
7 import "gonum.org/v1/gonum/blas"
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.
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.
18 // tau must have length at least n-1, and Dorgtr will panic otherwise.
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
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)
32 if len(work) < lwork {
35 if lwork < n-1 && lwork != -1 {
38 upper := uplo == blas.Upper
39 if !upper && uplo != blas.Lower {
50 nb = impl.Ilaenv(1, "DORGQL", " ", n-1, n-1, n-1, -1)
52 nb = impl.Ilaenv(1, "DORGQR", " ", n-1, n-1, n-1, -1)
54 lworkopt := max(1, n-1) * nb
56 work[0] = float64(lworkopt)
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
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]
71 for i := 0; i < n-1; i++ {
76 // Generate Q[0:n-1, 0:n-1].
77 impl.Dorgql(n-1, n-1, n-1, a, lda, tau, work, lwork)
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
83 for j := n - 1; j > 0; j-- {
85 for i := j + 1; i < n; i++ {
86 a[i*lda+j] = a[i*lda+j-1]
90 for i := 1; i < n; i++ {
94 // Generate Q[1:n, 1:n].
95 impl.Dorgqr(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork)
98 work[0] = float64(lworkopt)