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