github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dorghr.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  // Dorghr generates an n×n orthogonal matrix Q which is defined as the product
     8  // of ihi-ilo elementary reflectors:
     9  //  Q = H_{ilo} H_{ilo+1} ... H_{ihi-1}.
    10  //
    11  // a and lda represent an n×n matrix that contains the elementary reflectors, as
    12  // returned by Dgehrd. On return, a is overwritten by the n×n orthogonal matrix
    13  // Q. Q will be equal to the identity matrix except in the submatrix
    14  // Q[ilo+1:ihi+1,ilo+1:ihi+1].
    15  //
    16  // ilo and ihi must have the same values as in the previous call of Dgehrd. It
    17  // must hold that
    18  //  0 <= ilo <= ihi < n,  if n > 0,
    19  //  ilo = 0, ihi = -1,    if n == 0.
    20  //
    21  // tau contains the scalar factors of the elementary reflectors, as returned by
    22  // Dgehrd. tau must have length n-1.
    23  //
    24  // work must have length at least max(1,lwork) and lwork must be at least
    25  // ihi-ilo. For optimum performance lwork must be at least (ihi-ilo)*nb where nb
    26  // is the optimal blocksize. On return, work[0] will contain the optimal value
    27  // of lwork.
    28  //
    29  // If lwork == -1, instead of performing Dorghr, only the optimal value of lwork
    30  // will be stored into work[0].
    31  //
    32  // If any requirement on input sizes is not met, Dorghr will panic.
    33  //
    34  // Dorghr is an internal routine. It is exported for testing purposes.
    35  func (impl Implementation) Dorghr(n, ilo, ihi int, a []float64, lda int, tau, work []float64, lwork int) {
    36  	checkMatrix(n, n, a, lda)
    37  	nh := ihi - ilo
    38  	switch {
    39  	case ilo < 0 || max(1, n) <= ilo:
    40  		panic(badIlo)
    41  	case ihi < min(ilo, n-1) || n <= ihi:
    42  		panic(badIhi)
    43  	case lwork < max(1, nh) && lwork != -1:
    44  		panic(badWork)
    45  	case len(work) < max(1, lwork):
    46  		panic(shortWork)
    47  	}
    48  
    49  	lwkopt := max(1, nh) * impl.Ilaenv(1, "DORGQR", " ", nh, nh, nh, -1)
    50  	if lwork == -1 {
    51  		work[0] = float64(lwkopt)
    52  		return
    53  	}
    54  
    55  	// Quick return if possible.
    56  	if n == 0 {
    57  		work[0] = 1
    58  		return
    59  	}
    60  
    61  	// Shift the vectors which define the elementary reflectors one column
    62  	// to the right.
    63  	for i := ilo + 2; i < ihi+1; i++ {
    64  		copy(a[i*lda+ilo+1:i*lda+i], a[i*lda+ilo:i*lda+i-1])
    65  	}
    66  	// Set the first ilo+1 and the last n-ihi-1 rows and columns to those of
    67  	// the identity matrix.
    68  	for i := 0; i < ilo+1; i++ {
    69  		for j := 0; j < n; j++ {
    70  			a[i*lda+j] = 0
    71  		}
    72  		a[i*lda+i] = 1
    73  	}
    74  	for i := ilo + 1; i < ihi+1; i++ {
    75  		for j := 0; j <= ilo; j++ {
    76  			a[i*lda+j] = 0
    77  		}
    78  		for j := i; j < n; j++ {
    79  			a[i*lda+j] = 0
    80  		}
    81  	}
    82  	for i := ihi + 1; i < n; i++ {
    83  		for j := 0; j < n; j++ {
    84  			a[i*lda+j] = 0
    85  		}
    86  		a[i*lda+i] = 1
    87  	}
    88  	if nh > 0 {
    89  		// Generate Q[ilo+1:ihi+1,ilo+1:ihi+1].
    90  		impl.Dorgqr(nh, nh, nh, a[(ilo+1)*lda+ilo+1:], lda, tau[ilo:ihi], work, lwork)
    91  	}
    92  	work[0] = float64(lwkopt)
    93  }