github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dgehd2.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  // Dgehd2 reduces a block of a general n×n matrix A to upper Hessenberg form H
    10  // by an orthogonal similarity transformation Q^T * A * Q = H.
    11  //
    12  // The matrix Q is represented as a product of (ihi-ilo) elementary
    13  // reflectors
    14  //  Q = H_{ilo} H_{ilo+1} ... H_{ihi-1}.
    15  // Each H_i has the form
    16  //  H_i = I - tau[i] * v * v^T
    17  // where v is a real vector with v[0:i+1] = 0, v[i+1] = 1 and v[ihi+1:n] = 0.
    18  // v[i+2:ihi+1] is stored on exit in A[i+2:ihi+1,i].
    19  //
    20  // On entry, a contains the n×n general matrix to be reduced. On return, the
    21  // upper triangle and the first subdiagonal of A are overwritten with the upper
    22  // Hessenberg matrix H, and the elements below the first subdiagonal, with the
    23  // slice tau, represent the orthogonal matrix Q as a product of elementary
    24  // reflectors.
    25  //
    26  // The contents of A are illustrated by the following example, with n = 7, ilo =
    27  // 1 and ihi = 5.
    28  // On entry,
    29  //  [ a   a   a   a   a   a   a ]
    30  //  [     a   a   a   a   a   a ]
    31  //  [     a   a   a   a   a   a ]
    32  //  [     a   a   a   a   a   a ]
    33  //  [     a   a   a   a   a   a ]
    34  //  [     a   a   a   a   a   a ]
    35  //  [                         a ]
    36  // on return,
    37  //  [ a   a   h   h   h   h   a ]
    38  //  [     a   h   h   h   h   a ]
    39  //  [     h   h   h   h   h   h ]
    40  //  [     v1  h   h   h   h   h ]
    41  //  [     v1  v2  h   h   h   h ]
    42  //  [     v1  v2  v3  h   h   h ]
    43  //  [                         a ]
    44  // where a denotes an element of the original matrix A, h denotes a
    45  // modified element of the upper Hessenberg matrix H, and vi denotes an
    46  // element of the vector defining H_i.
    47  //
    48  // ilo and ihi determine the block of A that will be reduced to upper Hessenberg
    49  // form. It must hold that 0 <= ilo <= ihi <= max(0, n-1), otherwise Dgehd2 will
    50  // panic.
    51  //
    52  // On return, tau will contain the scalar factors of the elementary reflectors.
    53  // It must have length equal to n-1, otherwise Dgehd2 will panic.
    54  //
    55  // work must have length at least n, otherwise Dgehd2 will panic.
    56  //
    57  // Dgehd2 is an internal routine. It is exported for testing purposes.
    58  func (impl Implementation) Dgehd2(n, ilo, ihi int, a []float64, lda int, tau, work []float64) {
    59  	checkMatrix(n, n, a, lda)
    60  	switch {
    61  	case ilo < 0 || ilo > max(0, n-1):
    62  		panic(badIlo)
    63  	case ihi < min(ilo, n-1) || ihi >= n:
    64  		panic(badIhi)
    65  	case len(tau) != n-1:
    66  		panic(badTau)
    67  	case len(work) < n:
    68  		panic(badWork)
    69  	}
    70  
    71  	for i := ilo; i < ihi; i++ {
    72  		// Compute elementary reflector H_i to annihilate A[i+2:ihi+1,i].
    73  		var aii float64
    74  		aii, tau[i] = impl.Dlarfg(ihi-i, a[(i+1)*lda+i], a[min(i+2, n-1)*lda+i:], lda)
    75  		a[(i+1)*lda+i] = 1
    76  
    77  		// Apply H_i to A[0:ihi+1,i+1:ihi+1] from the right.
    78  		impl.Dlarf(blas.Right, ihi+1, ihi-i, a[(i+1)*lda+i:], lda, tau[i], a[i+1:], lda, work)
    79  
    80  		// Apply H_i to A[i+1:ihi+1,i+1:n] from the left.
    81  		impl.Dlarf(blas.Left, ihi-i, n-i-1, a[(i+1)*lda+i:], lda, tau[i], a[(i+1)*lda+i+1:], lda, work)
    82  		a[(i+1)*lda+i] = aii
    83  	}
    84  }