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.
8 "gonum.org/v1/gonum/blas"
9 "gonum.org/v1/gonum/lapack"
12 // Dorgql generates the m×n matrix Q with orthonormal columns defined as the
13 // last n columns of a product of k elementary reflectors of order m
14 // Q = H_{k-1} * ... * H_1 * H_0.
18 // and Dorgql will panic otherwise.
20 // On entry, the (n-k+i)-th column of A must contain the vector which defines
21 // the elementary reflector H_i, for i=0,...,k-1, and tau[i] must contain its
22 // scalar factor. On return, a contains the m×n matrix Q.
24 // tau must have length at least k, and Dorgql will panic otherwise.
26 // work must have length at least max(1,lwork), and lwork must be at least
27 // max(1,n), otherwise Dorgql will panic. For optimum performance lwork must
28 // be a sufficiently large multiple of n.
30 // If lwork == -1, instead of computing Dorgql the optimal work length is stored
33 // Dorgql is an internal routine. It is exported for testing purposes.
34 func (impl Implementation) Dorgql(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
44 case lwork < max(1, n) && lwork != -1:
46 case len(work) < lwork:
50 checkMatrix(m, n, a, lda)
61 nb := impl.Ilaenv(1, "DORGQL", " ", m, n, k, -1)
63 work[0] = float64(n * nb)
71 // Determine when to cross over from blocked to unblocked code.
72 nx = max(0, impl.Ilaenv(3, "DORGQL", " ", m, n, k, -1))
74 // Determine if workspace is large enough for blocked code.
77 // Not enough workspace to use optimal nb: reduce nb and determine
78 // the minimum value of nb.
80 nbmin = max(2, impl.Ilaenv(2, "DORGQL", " ", m, n, k, -1))
87 if nb >= nbmin && nb < k && nx < k {
88 // Use blocked code after the first block. The last kk columns are handled
89 // by the block method.
90 kk = min(k, ((k-nx+nb-1)/nb)*nb)
92 // Set A(m-kk:m, 0:n-kk) to zero.
93 for i := m - kk; i < m; i++ {
94 for j := 0; j < n-kk; j++ {
100 // Use unblocked code for the first or only block.
101 impl.Dorg2l(m-kk, n-kk, k-kk, a, lda, tau, work)
104 for i := k - kk; i < k; i += nb {
107 // Form the triangular factor of the block reflector
108 // H = H_{i+ib-1} * ... * H_{i+1} * H_i.
109 impl.Dlarft(lapack.Backward, lapack.ColumnWise, m-k+i+ib, ib,
110 a[n-k+i:], lda, tau[i:], work, ldwork)
112 // Apply H to A[0:m-k+i+ib, 0:n-k+i] from the left.
113 impl.Dlarfb(blas.Left, blas.NoTrans, lapack.Backward, lapack.ColumnWise,
114 m-k+i+ib, n-k+i, ib, a[n-k+i:], lda, work, ldwork,
115 a, lda, work[ib*ldwork:], ldwork)
118 // Apply H to rows 0:m-k+i+ib of current block.
119 impl.Dorg2l(m-k+i+ib, ib, ib, a[n-k+i:], lda, tau[i:], work)
121 // Set rows m-k+i+ib:m of current block to zero.
122 for j := n - k + i; j < n-k+i+ib; j++ {
123 for l := m - k + i + ib; l < m; l++ {
129 work[0] = float64(iws)