github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/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 gonum
     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  			safmn2 := math.Pow(dlamchB, math.Log(dlamchS/dlamchE)/math.Log(dlamchB)/2)
    70  			safmx2 := 1 / safmn2
    71  			sigma := b + c
    72  		loop:
    73  			for iter := 0; iter < 20; iter++ {
    74  				scale = math.Max(math.Abs(temp), math.Abs(sigma))
    75  				switch {
    76  				case scale >= safmx2:
    77  					sigma *= safmn2
    78  					temp *= safmn2
    79  				case scale <= safmn2:
    80  					sigma *= safmx2
    81  					temp *= safmx2
    82  				default:
    83  					break loop
    84  				}
    85  			}
    86  			p = temp / 2
    87  			tau := impl.Dlapy2(sigma, temp)
    88  			cs = math.Sqrt((1 + math.Abs(sigma)/tau) / 2)
    89  			sn = -p / (tau * cs)
    90  			if sigma < 0 {
    91  				sn *= -1
    92  			}
    93  			// Compute [ aa bb ] = [ a b ] [ cs -sn ]
    94  			//         [ cc dd ]   [ c d ] [ sn  cs ]
    95  			aa = a*cs + b*sn
    96  			bb = -a*sn + b*cs
    97  			cc = c*cs + d*sn
    98  			dd = -c*sn + d*cs
    99  			// Compute [ a b ] = [ cs sn ] [ aa bb ]
   100  			//         [ c d ]   [-sn cs ] [ cc dd ]
   101  			a = aa*cs + cc*sn
   102  			b = bb*cs + dd*sn
   103  			c = -aa*sn + cc*cs
   104  			d = -bb*sn + dd*cs
   105  
   106  			temp = (a + d) / 2
   107  			aa = temp
   108  			bb = b
   109  			cc = c
   110  			dd = temp
   111  
   112  			if cc != 0 {
   113  				if bb != 0 {
   114  					if math.Signbit(bb) == math.Signbit(cc) {
   115  						// Real eigenvalues, reduce to
   116  						// upper triangular form.
   117  						sab := math.Sqrt(math.Abs(bb))
   118  						sac := math.Sqrt(math.Abs(cc))
   119  						p = sab * sac
   120  						if cc < 0 {
   121  							p *= -1
   122  						}
   123  						tau = 1 / math.Sqrt(math.Abs(bb+cc))
   124  						aa = temp + p
   125  						bb = bb - cc
   126  						cc = 0
   127  						dd = temp - p
   128  						cs1 := sab * tau
   129  						sn1 := sac * tau
   130  						cs, sn = cs*cs1-sn*sn1, cs*sn1+sn*cs1
   131  					}
   132  				} else {
   133  					bb = -cc
   134  					cc = 0
   135  					cs, sn = -sn, cs
   136  				}
   137  			}
   138  		}
   139  	}
   140  
   141  	// Store eigenvalues in (rt1r,rt1i) and (rt2r,rt2i).
   142  	rt1r = aa
   143  	rt2r = dd
   144  	if cc != 0 {
   145  		rt1i = math.Sqrt(math.Abs(bb)) * math.Sqrt(math.Abs(cc))
   146  		rt2i = -rt1i
   147  	}
   148  	return
   149  }