github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dgelqf.go (about)

     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  	"github.com/jingcheng-WU/gonum/blas"
     9  	"github.com/jingcheng-WU/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  	switch {
    25  	case m < 0:
    26  		panic(mLT0)
    27  	case n < 0:
    28  		panic(nLT0)
    29  	case lda < max(1, n):
    30  		panic(badLdA)
    31  	case lwork < max(1, m) && lwork != -1:
    32  		panic(badLWork)
    33  	case len(work) < max(1, lwork):
    34  		panic(shortWork)
    35  	}
    36  
    37  	k := min(m, n)
    38  	if k == 0 {
    39  		work[0] = 1
    40  		return
    41  	}
    42  
    43  	nb := impl.Ilaenv(1, "DGELQF", " ", m, n, -1, -1)
    44  	if lwork == -1 {
    45  		work[0] = float64(m * nb)
    46  		return
    47  	}
    48  
    49  	if len(a) < (m-1)*lda+n {
    50  		panic(shortA)
    51  	}
    52  	if len(tau) < k {
    53  		panic(shortTau)
    54  	}
    55  
    56  	// Find the optimal blocking size based on the size of available memory
    57  	// and optimal machine parameters.
    58  	nbmin := 2
    59  	var nx int
    60  	iws := m
    61  	if 1 < nb && nb < k {
    62  		nx = max(0, impl.Ilaenv(3, "DGELQF", " ", m, n, -1, -1))
    63  		if nx < k {
    64  			iws = m * nb
    65  			if lwork < iws {
    66  				nb = lwork / m
    67  				nbmin = max(2, impl.Ilaenv(2, "DGELQF", " ", m, n, -1, -1))
    68  			}
    69  		}
    70  	}
    71  	ldwork := nb
    72  	// Computed blocked LQ factorization.
    73  	var i int
    74  	if nbmin <= nb && nb < k && nx < k {
    75  		for i = 0; i < k-nx; i += nb {
    76  			ib := min(k-i, nb)
    77  			impl.Dgelq2(ib, n-i, a[i*lda+i:], lda, tau[i:], work)
    78  			if i+ib < m {
    79  				impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
    80  					a[i*lda+i:], lda,
    81  					tau[i:],
    82  					work, ldwork)
    83  				impl.Dlarfb(blas.Right, blas.NoTrans, lapack.Forward, lapack.RowWise,
    84  					m-i-ib, n-i, ib,
    85  					a[i*lda+i:], lda,
    86  					work, ldwork,
    87  					a[(i+ib)*lda+i:], lda,
    88  					work[ib*ldwork:], ldwork)
    89  			}
    90  		}
    91  	}
    92  	// Perform unblocked LQ factorization on the remainder.
    93  	if i < k {
    94  		impl.Dgelq2(m-i, n-i, a[i*lda+i:], lda, tau[i:], work)
    95  	}
    96  	work[0] = float64(iws)
    97  }