github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dorgtr.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 "github.com/gonum/blas"
     8  
     9  // Dorgtr generates a real orthogonal matrix Q which is defined as the product
    10  // of n-1 elementary reflectors of order n as returned by Dsytrd.
    11  //
    12  // The construction of Q depends on the value of uplo:
    13  //  Q = H_{n-1} * ... * H_1 * H_0  if uplo == blas.Upper
    14  //  Q = H_0 * H_1 * ... * H_{n-1}  if uplo == blas.Lower
    15  // where H_i is constructed from the elementary reflectors as computed by Dsytrd.
    16  // See the documentation for Dsytrd for more information.
    17  //
    18  // tau must have length at least n-1, and Dorgtr will panic otherwise.
    19  //
    20  // work is temporary storage, and lwork specifies the usable memory length. At
    21  // minimum, lwork >= max(1,n-1), and Dorgtr will panic otherwise. The amount of blocking
    22  // is limited by the usable length.
    23  // If lwork == -1, instead of computing Dorgtr the optimal work length is stored
    24  // into work[0].
    25  //
    26  // Dorgtr is an internal routine. It is exported for testing purposes.
    27  func (impl Implementation) Dorgtr(uplo blas.Uplo, n int, a []float64, lda int, tau, work []float64, lwork int) {
    28  	checkMatrix(n, n, a, lda)
    29  	if len(tau) < n-1 {
    30  		panic(badTau)
    31  	}
    32  	if len(work) < lwork {
    33  		panic(badWork)
    34  	}
    35  	if lwork < n-1 && lwork != -1 {
    36  		panic(badWork)
    37  	}
    38  	upper := uplo == blas.Upper
    39  	if !upper && uplo != blas.Lower {
    40  		panic(badUplo)
    41  	}
    42  
    43  	if n == 0 {
    44  		work[0] = 1
    45  		return
    46  	}
    47  
    48  	var nb int
    49  	if upper {
    50  		nb = impl.Ilaenv(1, "DORGQL", " ", n-1, n-1, n-1, -1)
    51  	} else {
    52  		nb = impl.Ilaenv(1, "DORGQR", " ", n-1, n-1, n-1, -1)
    53  	}
    54  	lworkopt := max(1, n-1) * nb
    55  	if lwork == -1 {
    56  		work[0] = float64(lworkopt)
    57  		return
    58  	}
    59  
    60  	if upper {
    61  		// Q was determined by a call to Dsytrd with uplo == blas.Upper.
    62  		// Shift the vectors which define the elementary reflectors one column
    63  		// to the left, and set the last row and column of Q to those of the unit
    64  		// matrix.
    65  		for j := 0; j < n-1; j++ {
    66  			for i := 0; i < j; i++ {
    67  				a[i*lda+j] = a[i*lda+j+1]
    68  			}
    69  			a[(n-1)*lda+j] = 0
    70  		}
    71  		for i := 0; i < n-1; i++ {
    72  			a[i*lda+n-1] = 0
    73  		}
    74  		a[(n-1)*lda+n-1] = 1
    75  
    76  		// Generate Q[0:n-1, 0:n-1].
    77  		impl.Dorgql(n-1, n-1, n-1, a, lda, tau, work, lwork)
    78  	} else {
    79  		// Q was determined by a call to Dsytrd with uplo == blas.Upper.
    80  		// Shift the vectors which define the elementary reflectors one column
    81  		// to the right, and set the first row and column of Q to those of the unit
    82  		// matrix.
    83  		for j := n - 1; j > 0; j-- {
    84  			a[j] = 0
    85  			for i := j + 1; i < n; i++ {
    86  				a[i*lda+j] = a[i*lda+j-1]
    87  			}
    88  		}
    89  		a[0] = 1
    90  		for i := 1; i < n; i++ {
    91  			a[i*lda] = 0
    92  		}
    93  		if n > 1 {
    94  			// Generate Q[1:n, 1:n].
    95  			impl.Dorgqr(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork)
    96  		}
    97  	}
    98  	work[0] = float64(lworkopt)
    99  }