github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dlanv2.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 "math"
     8  
     9  // Dlanv2 computes the Schur factorization of a real 2×2 matrix:
    10  //  [ a b ] = [ cs -sn ] * [ aa bb ] * [ cs sn ]
    11  //  [ c d ]   [ sn  cs ]   [ cc dd ] * [-sn cs ]
    12  // If cc is zero, aa and dd are real eigenvalues of the matrix. Otherwise it
    13  // holds that aa = dd and bb*cc < 0, and aa ± sqrt(bb*cc) are complex conjugate
    14  // eigenvalues. The real and imaginary parts of the eigenvalues are returned in
    15  // (rt1r,rt1i) and (rt2r,rt2i).
    16  func (impl Implementation) Dlanv2(a, b, c, d float64) (aa, bb, cc, dd float64, rt1r, rt1i, rt2r, rt2i float64, cs, sn float64) {
    17  	switch {
    18  	case c == 0: // Matrix is already upper triangular.
    19  		aa = a
    20  		bb = b
    21  		cc = 0
    22  		dd = d
    23  		cs = 1
    24  		sn = 0
    25  	case b == 0: // Matrix is lower triangular, swap rows and columns.
    26  		aa = d
    27  		bb = -c
    28  		cc = 0
    29  		dd = a
    30  		cs = 0
    31  		sn = 1
    32  	case a == d && math.Signbit(b) != math.Signbit(c): // Matrix is already in the standard Schur form.
    33  		aa = a
    34  		bb = b
    35  		cc = c
    36  		dd = d
    37  		cs = 1
    38  		sn = 0
    39  	default:
    40  		temp := a - d
    41  		p := temp / 2
    42  		bcmax := math.Max(math.Abs(b), math.Abs(c))
    43  		bcmis := math.Min(math.Abs(b), math.Abs(c))
    44  		if b*c < 0 {
    45  			bcmis *= -1
    46  		}
    47  		scale := math.Max(math.Abs(p), bcmax)
    48  		z := p/scale*p + bcmax/scale*bcmis
    49  		eps := dlamchP
    50  
    51  		if z >= 4*eps {
    52  			// Real eigenvalues. Compute aa and dd.
    53  			if p > 0 {
    54  				z = p + math.Sqrt(scale)*math.Sqrt(z)
    55  			} else {
    56  				z = p - math.Sqrt(scale)*math.Sqrt(z)
    57  			}
    58  			aa = d + z
    59  			dd = d - bcmax/z*bcmis
    60  			// Compute bb and the rotation matrix.
    61  			tau := impl.Dlapy2(c, z)
    62  			cs = z / tau
    63  			sn = c / tau
    64  			bb = b - c
    65  			cc = 0
    66  		} else {
    67  			// Complex eigenvalues, or real (almost) equal eigenvalues.
    68  			// Make diagonal elements equal.
    69  			sigma := b + c
    70  			tau := impl.Dlapy2(sigma, temp)
    71  			cs = math.Sqrt((1 + math.Abs(sigma)/tau) / 2)
    72  			sn = -p / (tau * cs)
    73  			if sigma < 0 {
    74  				sn *= -1
    75  			}
    76  			// Compute [ aa bb ] = [ a b ] [ cs -sn ]
    77  			//         [ cc dd ]   [ c d ] [ sn  cs ]
    78  			aa = a*cs + b*sn
    79  			bb = -a*sn + b*cs
    80  			cc = c*cs + d*sn
    81  			dd = -c*sn + d*cs
    82  			// Compute [ a b ] = [ cs sn ] [ aa bb ]
    83  			//         [ c d ]   [-sn cs ] [ cc dd ]
    84  			a = aa*cs + cc*sn
    85  			b = bb*cs + dd*sn
    86  			c = -aa*sn + cc*cs
    87  			d = -bb*sn + dd*cs
    88  
    89  			temp = (a + d) / 2
    90  			aa = temp
    91  			bb = b
    92  			cc = c
    93  			dd = temp
    94  
    95  			if cc != 0 {
    96  				if bb != 0 {
    97  					if math.Signbit(bb) == math.Signbit(cc) {
    98  						// Real eigenvalues, reduce to
    99  						// upper triangular form.
   100  						sab := math.Sqrt(math.Abs(bb))
   101  						sac := math.Sqrt(math.Abs(cc))
   102  						p = sab * sac
   103  						if cc < 0 {
   104  							p *= -1
   105  						}
   106  						tau = 1 / math.Sqrt(math.Abs(bb+cc))
   107  						aa = temp + p
   108  						bb = bb - cc
   109  						cc = 0
   110  						dd = temp - p
   111  						cs1 := sab * tau
   112  						sn1 := sac * tau
   113  						cs, sn = cs*cs1-sn*sn1, cs*sn1+sn+cs1
   114  					}
   115  				} else {
   116  					bb = -cc
   117  					cc = 0
   118  					cs, sn = -sn, cs
   119  				}
   120  			}
   121  		}
   122  	}
   123  
   124  	// Store eigenvalues in (rt1r,rt1i) and (rt2r,rt2i).
   125  	rt1r = aa
   126  	rt2r = dd
   127  	if cc != 0 {
   128  		rt1i = math.Sqrt(math.Abs(bb)) * math.Sqrt(math.Abs(cc))
   129  		rt2i = -rt1i
   130  	}
   131  	return
   132  }