gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/gonum/dgghrd.go (about)

     1  // Copyright ©2023 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  	"gonum.org/v1/gonum/blas"
     9  	"gonum.org/v1/gonum/blas/blas64"
    10  	"gonum.org/v1/gonum/lapack"
    11  )
    12  
    13  // Dgghrd reduces a pair of real matrices (A,B) to generalized upper Hessenberg
    14  // form using orthogonal transformations, where A is a general matrix and B is
    15  // upper triangular.
    16  //
    17  // This subroutine simultaneously reduces A to a Hessenberg matrix H
    18  //
    19  //	Qᵀ*A*Z = H,
    20  //
    21  // and transforms B to another upper triangular matrix T
    22  //
    23  //	Qᵀ*B*Z = T.
    24  //
    25  // The orthogonal matrices Q and Z are determined as products of Givens
    26  // rotations. They may either be formed explicitly (lapack.OrthoExplicit), or
    27  // they may be postmultiplied into input matrices Q1 and Z1
    28  // (lapack.OrthoPostmul), so that
    29  //
    30  //	Q1 * A * Z1ᵀ = (Q1*Q) * H * (Z1*Z)ᵀ,
    31  //	Q1 * B * Z1ᵀ = (Q1*Q) * T * (Z1*Z)ᵀ.
    32  //
    33  // ilo and ihi determine the block of A that will be reduced. It must hold that
    34  //
    35  //   - 0 <= ilo <= ihi < n      if n > 0,
    36  //   - ilo == 0 and ihi == -1   if n == 0,
    37  //
    38  // otherwise Dgghrd will panic.
    39  //
    40  // Dgghrd is an internal routine. It is exported for testing purposes.
    41  func (impl Implementation) Dgghrd(compq, compz lapack.OrthoComp, n, ilo, ihi int, a []float64, lda int, b []float64, ldb int, q []float64, ldq int, z []float64, ldz int) {
    42  	switch {
    43  	case compq != lapack.OrthoNone && compq != lapack.OrthoExplicit && compq != lapack.OrthoPostmul:
    44  		panic(badOrthoComp)
    45  	case compz != lapack.OrthoNone && compz != lapack.OrthoExplicit && compz != lapack.OrthoPostmul:
    46  		panic(badOrthoComp)
    47  	case n < 0:
    48  		panic(nLT0)
    49  	case ilo < 0 || max(0, n-1) < ilo:
    50  		panic(badIlo)
    51  	case ihi < min(ilo, n-1) || n <= ihi:
    52  		panic(badIhi)
    53  	case lda < max(1, n):
    54  		panic(badLdA)
    55  	case ldb < max(1, n):
    56  		panic(badLdB)
    57  	case (compq != lapack.OrthoNone && ldq < n) || ldq < 1:
    58  		panic(badLdQ)
    59  	case (compz != lapack.OrthoNone && ldz < n) || ldz < 1:
    60  		panic(badLdZ)
    61  	}
    62  
    63  	// Quick return if possible.
    64  	if n == 0 {
    65  		return
    66  	}
    67  
    68  	switch {
    69  	case len(a) < (n-1)*lda+n:
    70  		panic(shortA)
    71  	case len(b) < (n-1)*ldb+n:
    72  		panic(shortB)
    73  	case compq != lapack.OrthoNone && len(q) < (n-1)*ldq+n:
    74  		panic(shortQ)
    75  	case compz != lapack.OrthoNone && len(z) < (n-1)*ldz+n:
    76  		panic(shortZ)
    77  	}
    78  
    79  	if compq == lapack.OrthoExplicit {
    80  		impl.Dlaset(blas.All, n, n, 0, 1, q, ldq)
    81  	}
    82  	if compz == lapack.OrthoExplicit {
    83  		impl.Dlaset(blas.All, n, n, 0, 1, z, ldz)
    84  	}
    85  
    86  	// Quick return if possible.
    87  	if n == 1 {
    88  		return
    89  	}
    90  
    91  	// Zero out lower triangle of B.
    92  	for i := 1; i < n; i++ {
    93  		for j := 0; j < i; j++ {
    94  			b[i*ldb+j] = 0
    95  		}
    96  	}
    97  	bi := blas64.Implementation()
    98  	// Reduce A and B.
    99  	for jcol := ilo; jcol <= ihi-2; jcol++ {
   100  		for jrow := ihi; jrow >= jcol+2; jrow-- {
   101  			// Step 1: rotate rows jrow-1, jrow to kill A[jrow,jcol].
   102  			var c, s float64
   103  			c, s, a[(jrow-1)*lda+jcol] = impl.Dlartg(a[(jrow-1)*lda+jcol], a[jrow*lda+jcol])
   104  			a[jrow*lda+jcol] = 0
   105  
   106  			bi.Drot(n-jcol-1, a[(jrow-1)*lda+jcol+1:], 1, a[jrow*lda+jcol+1:], 1, c, s)
   107  			bi.Drot(n+2-jrow-1, b[(jrow-1)*ldb+jrow-1:], 1, b[jrow*ldb+jrow-1:], 1, c, s)
   108  
   109  			if compq != lapack.OrthoNone {
   110  				bi.Drot(n, q[jrow-1:], ldq, q[jrow:], ldq, c, s)
   111  			}
   112  
   113  			// Step 2: rotate columns jrow, jrow-1 to kill B[jrow,jrow-1].
   114  			c, s, b[jrow*ldb+jrow] = impl.Dlartg(b[jrow*ldb+jrow], b[jrow*ldb+jrow-1])
   115  			b[jrow*ldb+jrow-1] = 0
   116  
   117  			bi.Drot(ihi+1, a[jrow:], lda, a[jrow-1:], lda, c, s)
   118  			bi.Drot(jrow, b[jrow:], ldb, b[jrow-1:], ldb, c, s)
   119  
   120  			if compz != lapack.OrthoNone {
   121  				bi.Drot(n, z[jrow:], ldz, z[jrow-1:], ldz, c, s)
   122  			}
   123  		}
   124  	}
   125  }