github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/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 gonum
     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  	nh := ihi - ilo
    37  	switch {
    38  	case ilo < 0 || max(1, n) <= ilo:
    39  		panic(badIlo)
    40  	case ihi < min(ilo, n-1) || n <= ihi:
    41  		panic(badIhi)
    42  	case lda < max(1, n):
    43  		panic(badLdA)
    44  	case lwork < max(1, nh) && lwork != -1:
    45  		panic(badLWork)
    46  	case len(work) < max(1, lwork):
    47  		panic(shortWork)
    48  	}
    49  
    50  	// Quick return if possible.
    51  	if n == 0 {
    52  		work[0] = 1
    53  		return
    54  	}
    55  
    56  	lwkopt := max(1, nh) * impl.Ilaenv(1, "DORGQR", " ", nh, nh, nh, -1)
    57  	if lwork == -1 {
    58  		work[0] = float64(lwkopt)
    59  		return
    60  	}
    61  
    62  	switch {
    63  	case len(a) < (n-1)*lda+n:
    64  		panic(shortA)
    65  	case len(tau) < n-1:
    66  		panic(shortTau)
    67  	}
    68  
    69  	// Shift the vectors which define the elementary reflectors one column
    70  	// to the right.
    71  	for i := ilo + 2; i < ihi+1; i++ {
    72  		copy(a[i*lda+ilo+1:i*lda+i], a[i*lda+ilo:i*lda+i-1])
    73  	}
    74  	// Set the first ilo+1 and the last n-ihi-1 rows and columns to those of
    75  	// the identity matrix.
    76  	for i := 0; i < ilo+1; i++ {
    77  		for j := 0; j < n; j++ {
    78  			a[i*lda+j] = 0
    79  		}
    80  		a[i*lda+i] = 1
    81  	}
    82  	for i := ilo + 1; i < ihi+1; i++ {
    83  		for j := 0; j <= ilo; j++ {
    84  			a[i*lda+j] = 0
    85  		}
    86  		for j := i; j < n; j++ {
    87  			a[i*lda+j] = 0
    88  		}
    89  	}
    90  	for i := ihi + 1; i < n; i++ {
    91  		for j := 0; j < n; j++ {
    92  			a[i*lda+j] = 0
    93  		}
    94  		a[i*lda+i] = 1
    95  	}
    96  	if nh > 0 {
    97  		// Generate Q[ilo+1:ihi+1,ilo+1:ihi+1].
    98  		impl.Dorgqr(nh, nh, nh, a[(ilo+1)*lda+ilo+1:], lda, tau[ilo:ihi], work, lwork)
    99  	}
   100  	work[0] = float64(lwkopt)
   101  }