gonum.org/v1/gonum@v0.14.0/lapack/gonum/dlaqr1.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  import "math"
     8  
     9  // Dlaqr1 sets v to a scalar multiple of the first column of the product
    10  //
    11  //	(H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)
    12  //
    13  // where H is a 2×2 or 3×3 matrix, I is the identity matrix of the same size,
    14  // and i is the imaginary unit. Scaling is done to avoid overflows and most
    15  // underflows.
    16  //
    17  // n is the order of H and must be either 2 or 3. It must hold that either sr1 =
    18  // sr2 and si1 = -si2, or si1 = si2 = 0. The length of v must be equal to n. If
    19  // any of these conditions is not met, Dlaqr1 will panic.
    20  //
    21  // Dlaqr1 is an internal routine. It is exported for testing purposes.
    22  func (impl Implementation) Dlaqr1(n int, h []float64, ldh int, sr1, si1, sr2, si2 float64, v []float64) {
    23  	switch {
    24  	case n != 2 && n != 3:
    25  		panic("lapack: n must be 2 or 3")
    26  	case ldh < n:
    27  		panic(badLdH)
    28  	case len(h) < (n-1)*ldh+n:
    29  		panic(shortH)
    30  	case !((sr1 == sr2 && si1 == -si2) || (si1 == 0 && si2 == 0)):
    31  		panic(badShifts)
    32  	case len(v) != n:
    33  		panic(shortV)
    34  	}
    35  
    36  	if n == 2 {
    37  		s := math.Abs(h[0]-sr2) + math.Abs(si2) + math.Abs(h[ldh])
    38  		if s == 0 {
    39  			v[0] = 0
    40  			v[1] = 0
    41  		} else {
    42  			h21s := h[ldh] / s
    43  			v[0] = h21s*h[1] + (h[0]-sr1)*((h[0]-sr2)/s) - si1*(si2/s)
    44  			v[1] = h21s * (h[0] + h[ldh+1] - sr1 - sr2)
    45  		}
    46  		return
    47  	}
    48  
    49  	s := math.Abs(h[0]-sr2) + math.Abs(si2) + math.Abs(h[ldh]) + math.Abs(h[2*ldh])
    50  	if s == 0 {
    51  		v[0] = 0
    52  		v[1] = 0
    53  		v[2] = 0
    54  	} else {
    55  		h21s := h[ldh] / s
    56  		h31s := h[2*ldh] / s
    57  		v[0] = (h[0]-sr1)*((h[0]-sr2)/s) - si1*(si2/s) + h[1]*h21s + h[2]*h31s
    58  		v[1] = h21s*(h[0]+h[ldh+1]-sr1-sr2) + h[ldh+2]*h31s
    59  		v[2] = h31s*(h[0]+h[2*ldh+2]-sr1-sr2) + h21s*h[2*ldh+1]
    60  	}
    61  }