gonum.org/v1/gonum@v0.14.0/lapack/gonum/dorgr2.go (about)

     1  // Copyright ©2021 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/blas/blas64"
    10  )
    11  
    12  // Dorgr2 generates an m×n real matrix Q with orthonormal rows, which is defined
    13  // as the last m rows of a product of k elementary reflectors of order n
    14  //
    15  //	Q = H_0 * H_1 * ... * H_{k-1}
    16  //
    17  // as returned by Dgerqf.
    18  //
    19  // On entry, the (m-k+i)-th row of A must contain the vector which defines the
    20  // elementary reflector H_i, for i = 0,1,...,k, as returned by Dgerqf. On
    21  // return, A will contain the m×n matrix Q.
    22  //
    23  // The i-th element of tau must contain the scalar factor of the elementary
    24  // reflector H_i, as returned by Dgerqf.
    25  //
    26  // It must hold that
    27  //
    28  //	n >= m >= k >= 0,
    29  //
    30  // the length of tau must be k and the length of work must be m, otherwise
    31  // Dorgr2 will panic.
    32  //
    33  // Dorgr2 is an internal routine. It is exported for testing purposes.
    34  func (impl Implementation) Dorgr2(m, n, k int, a []float64, lda int, tau, work []float64) {
    35  	switch {
    36  	case k < 0:
    37  		panic(kLT0)
    38  	case m < k:
    39  		panic(kGTM)
    40  	case n < m:
    41  		panic(mGTN)
    42  	case lda < max(1, n):
    43  		panic(badLdA)
    44  	}
    45  
    46  	// Quick return if possible.
    47  	if m == 0 {
    48  		return
    49  	}
    50  
    51  	switch {
    52  	case len(tau) != k:
    53  		panic(badLenTau)
    54  	case len(a) < (m-1)*lda+n:
    55  		panic(shortA)
    56  	case len(work) < m:
    57  		panic(shortWork)
    58  	}
    59  
    60  	// Initialise rows 0:m-k to rows of the unit matrix.
    61  	for l := 0; l < m-k; l++ {
    62  		row := a[l*lda : l*lda+n]
    63  		for j := range row {
    64  			row[j] = 0
    65  		}
    66  		a[l*lda+n-m+l] = 1
    67  	}
    68  	bi := blas64.Implementation()
    69  	for i := 0; i < k; i++ {
    70  		ii := m - k + i
    71  
    72  		// Apply H_i to A[0:m-k+i+1, 0:n-k+i+1] from the right.
    73  		a[ii*lda+n-m+ii] = 1
    74  		impl.Dlarf(blas.Right, ii, n-m+ii+1, a[ii*lda:], 1, tau[i], a, lda, work)
    75  		bi.Dscal(n-m+ii, -tau[i], a[ii*lda:], 1)
    76  		a[ii*lda+n-m+ii] = 1 - tau[i]
    77  
    78  		// Set A[m-k+i, n-k+i:n] to zero.
    79  		for l := n - m + ii + 1; l < n; l++ {
    80  			a[ii*lda+l] = 0
    81  		}
    82  	}
    83  }