github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"github.com/gonum/blas"
     9  	"github.com/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  }