vendor/gonum.org/v1/gonum/lapack/gonum/dgelqf.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 (
   8         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/lapack"
  10 )
  11
  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.
  15 //
  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].
  21 //
  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)
  26         if lwork == -1 {
  27                 work[0] = float64(lworkopt)
  28                 return
  29         }
  30         checkMatrix(m, n, a, lda)
  31         if len(work) < lwork {
  32                 panic(shortWork)
  33         }
  34         if lwork < m {
  35                 panic(badWork)
  36         }
  37         k := min(m, n)
  38         if len(tau) < k {
  39                 panic(badTau)
  40         }
  41         if k == 0 {
  42                 return
  43         }
  44         // Find the optimal blocking size based on the size of available memory
  45         // and optimal machine parameters.
  46         nbmin := 2
  47         var nx int
  48         iws := m
  49         ldwork := nb
  50         if nb > 1 && k > nb {
  51                 nx = max(0, impl.Ilaenv(3, "DGELQF", " ", m, n, -1, -1))
  52                 if nx < k {
  53                         iws = m * nb
  54                         if lwork < iws {
  55                                 nb = lwork / m
  56                                 nbmin = max(2, impl.Ilaenv(2, "DGELQF", " ", m, n, -1, -1))
  57                         }
  58                 }
  59         }
  60         // Computed blocked LQ factorization.
  61         var i int
  62         if nb >= nbmin && nb < k && nx < k {
  63                 for i = 0; i < k-nx; i += nb {
  64                         ib := min(k-i, nb)
  65                         impl.Dgelq2(ib, n-i, a[i*lda+i:], lda, tau[i:], work)
  66                         if i+ib < m {
  67                                 impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
  68                                         a[i*lda+i:], lda,
  69                                         tau[i:],
  70                                         work, ldwork)
  71                                 impl.Dlarfb(blas.Right, blas.NoTrans, lapack.Forward, lapack.RowWise,
  72                                         m-i-ib, n-i, ib,
  73                                         a[i*lda+i:], lda,
  74                                         work, ldwork,
  75                                         a[(i+ib)*lda+i:], lda,
  76                                         work[ib*ldwork:], ldwork)
  77                         }
  78                 }
  79         }
  80         // Perform unblocked LQ factorization on the remainder.
  81         if i < k {
  82                 impl.Dgelq2(m-i, n-i, a[i*lda+i:], lda, tau[i:], work)
  83         }
  84 }