#! /usr/bin/env python # def drotg ( a, b ): #*****************************************************************************80 # ## DROTG constructs a Givens plane rotation. # # Discussion: # # Given values A and B, this routine computes # # SIGMA = sign ( A ) if abs ( A ) > abs ( B ) # = sign ( B ) if abs ( A ) <= abs ( B ) # # R = SIGMA * ( A * A + B * B ) # # C = A / R if R is not 0 # = 1 if R is 0 # # S = B / R if R is not 0, # 0 if R is 0. # # The computed numbers then satisfy the equation # # ( C S ) ( A ) = ( R ) # ( -S C ) ( B ) = ( 0 ) # # The routine also computes # # Z = S if abs ( A ) > abs ( B ), # = 1 / C if abs ( A ) <= abs ( B ) and C is not 0, # = 1 if C is 0. # # The single value Z encodes C and S, and hence the rotation: # # If Z = 1, set C = 0 and S = 1 # If abs ( Z ) < 1, set C = sqrt ( 1 - Z * Z ) and S = Z # if abs ( Z ) > 1, set C = 1/ Z and S = sqrt ( 1 - C * C ) # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 06 September 2016 # # Author: # # Python version by John Burkardt. # # Reference: # # Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, # LINPACK User's Guide, # SIAM, 1979. # # Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, # Basic Linear Algebra Subprograms for Fortran Usage, # Algorithm 539, # ACM Transactions on Mathematical Software, # Volume 5, Number 3, September 1979, pages 308-323. # # Parameters: # # Input, real A, B, the values A and B. # # Output, real C, S, the cosine and sine of the Givens rotation. # # Output, real R, Z, the values R and Z. # import numpy as np if ( abs ( b ) < abs ( a ) ): roe = a else: roe = b scale = abs ( a ) + abs ( b ) if ( scale == 0.0 ): c = 1.0 s = 0.0 r = 0.0 else: r = scale * np.sqrt ( ( a / scale ) ** 2 + ( b / scale ) ** 2 ) if ( roe < 0.0 ): r = - r c = a / r s = b / r if ( 0.0 < abs ( c ) and abs ( c ) <= s ): z = 1.0 / c else: z = s return c, s, r, z def drotg_test ( ): #*****************************************************************************80 # ## DROTG_TEST tests DROTG. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 06 September 2016 # # Author: # # John Burkardt # import platform test_num = 5 print ( '' ) print ( 'DROTG_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' DROTG generates a real Givens rotation' ) print ( ' ( C S ) * ( A ) = ( R )' ) print ( ' ( -S C ) ( B ) ( 0 )' ) print ( '' ) seed = 123456789 for test in range ( 0, test_num ): a, seed = r8_uniform_01 ( seed ) b, seed = r8_uniform_01 ( seed ) c, s, r, z = drotg ( a, b ) print ( '' ) print ( ' A = %g B = %g' % ( a, b ) ) print ( ' C = %g S = %g' % ( c, s ) ) print ( ' R = %g Z = %g' % ( r, z ) ) print ( ' C*A+S*B = %g' % ( c * a + s * b ) ) print ( ' -S*A+C*B = %g' % ( -s * a + c * b ) ) # # Terminate. # print ( '' ) print ( 'DROTG_TEST' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) drotg_test ( ) timestamp ( )