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