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.
8 "gonum.org/v1/gonum/blas"
9 "gonum.org/v1/gonum/lapack"
12 // Dgelqf computes the LQ factorization of the m×n matrix A using a blocked
13 // algorithm. See the documentation for Dgelq2 for a description of the
14 // parameters at entry and exit.
16 // work is temporary storage, and lwork specifies the usable memory length.
17 // At minimum, lwork >= m, and this function will panic otherwise.
18 // Dgelqf is a blocked LQ factorization, but the block size is limited
19 // by the temporary space available. If lwork == -1, instead of performing Dgelqf,
20 // the optimal work length will be stored into work[0].
22 // tau must have length at least min(m,n), and this function will panic otherwise.
23 func (impl Implementation) Dgelqf(m, n int, a []float64, lda int, tau, work []float64, lwork int) {
24 nb := impl.Ilaenv(1, "DGELQF", " ", m, n, -1, -1)
25 lworkopt := m * max(nb, 1)
27 work[0] = float64(lworkopt)
30 checkMatrix(m, n, a, lda)
31 if len(work) < lwork {
44 // Find the optimal blocking size based on the size of available memory
45 // and optimal machine parameters.
51 nx = max(0, impl.Ilaenv(3, "DGELQF", " ", m, n, -1, -1))
56 nbmin = max(2, impl.Ilaenv(2, "DGELQF", " ", m, n, -1, -1))
60 // Computed blocked LQ factorization.
62 if nb >= nbmin && nb < k && nx < k {
63 for i = 0; i < k-nx; i += nb {
65 impl.Dgelq2(ib, n-i, a[i*lda+i:], lda, tau[i:], work)
67 impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
71 impl.Dlarfb(blas.Right, blas.NoTrans, lapack.Forward, lapack.RowWise,
75 a[(i+ib)*lda+i:], lda,
76 work[ib*ldwork:], ldwork)
80 // Perform unblocked LQ factorization on the remainder.
82 impl.Dgelq2(m-i, n-i, a[i*lda+i:], lda, tau[i:], work)