github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dorg2l.go (about)

     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.
     4  
     5  package native
     6  
     7  import (
     8  	"github.com/gonum/blas"
     9  	"github.com/gonum/blas/blas64"
    10  )
    11  
    12  // Dorg2l generates an m×n matrix Q with orthonormal columns which is defined
    13  // as the last n columns of a product of k elementary reflectors of order m.
    14  //  Q = H_{k-1} * ... * H_1 * H_0
    15  // See Dgelqf for more information. It must be that m >= n >= k.
    16  //
    17  // tau contains the scalar reflectors computed by Dgeqlf. tau must have length
    18  // at least k, and Dorg2l will panic otherwise.
    19  //
    20  // work contains temporary memory, and must have length at least n. Dorg2l will
    21  // panic otherwise.
    22  //
    23  // Dorg2l is an internal routine. It is exported for testing purposes.
    24  func (impl Implementation) Dorg2l(m, n, k int, a []float64, lda int, tau, work []float64) {
    25  	checkMatrix(m, n, a, lda)
    26  	if len(tau) < k {
    27  		panic(badTau)
    28  	}
    29  	if len(work) < n {
    30  		panic(badWork)
    31  	}
    32  	if m < n {
    33  		panic(mLTN)
    34  	}
    35  	if k > n {
    36  		panic(kGTN)
    37  	}
    38  	if n == 0 {
    39  		return
    40  	}
    41  
    42  	// Initialize columns 0:n-k to columns of the unit matrix.
    43  	for j := 0; j < n-k; j++ {
    44  		for l := 0; l < m; l++ {
    45  			a[l*lda+j] = 0
    46  		}
    47  		a[(m-n+j)*lda+j] = 1
    48  	}
    49  
    50  	bi := blas64.Implementation()
    51  	for i := 0; i < k; i++ {
    52  		ii := n - k + i
    53  
    54  		// Apply H_i to A[0:m-k+i, 0:n-k+i] from the left.
    55  		a[(m-n+ii)*lda+ii] = 1
    56  		impl.Dlarf(blas.Left, m-n+ii+1, ii, a[ii:], lda, tau[i], a, lda, work)
    57  		bi.Dscal(m-n+ii, -tau[i], a[ii:], lda)
    58  		a[(m-n+ii)*lda+ii] = 1 - tau[i]
    59  
    60  		// Set A[m-k+i:m, n-k+i+1] to zero.
    61  		for l := m - n + ii + 1; l < m; l++ {
    62  			a[l*lda+ii] = 0
    63  		}
    64  	}
    65  }