#! /usr/bin/env python # def drot ( n, x, incx, y, incy, c, s ): #*****************************************************************************80 # ## DROT applies a plane rotation. # # Discussion: # # Note that the FORTRAN version of this function allowed users to pass in # X and Y data that was noncontiguous, (such as rows of a FORTRAN matrix). # The rotated data overwrote the input data, and so it might therefore # also be noncontiguous. # # This function does not assume that the output overwrites the input, # and treats the output vectors as new items of length exactly N. # # Note, moreover, that Python does NOT allow a matrix to be treated as a # vector in quite the simple way that FORTRAN does. # # 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, integer N, the number of entries in the vectors. # # Input, real X(INCX*N), one of the vectors to be rotated. # # Input, integer INCX, the increment between successive entries of X. # # Input, real Y(INCX*N), one of the vectors to be rotated. # # Input, integer INCY, the increment between successive elements of Y. # # Input, real C, S, parameters (presumably the cosine and # sine of some angle) that define a plane rotation. # # Output, real XR(N), the rotated vector. # # Output, real YR(N), the rotated vector. # import numpy as np xr = np.zeros ( n ) yr = np.zeros ( n ) if ( n <= 0 ): pass elif ( incx == 1 and incy == 1 ): xr[0:n] = c * x[0:n] + s * y[0:n] yr[0:n] = c * y[0:n] - s * x[0:n] else: if ( 0 <= incx ): ix = 0 else: ix = ( - n + 1 ) * incx if ( 0 <= incy ): iy = 0 else: iy = ( - n + 1 ) * incy for i in range ( 0, n ): xr[ix] = c * x[ix] + s * y[iy] yr[iy] = c * y[iy] - s * x[ix] ix = ix + incx iy = iy + incy return xr, yr def drot_test ( ): #*****************************************************************************80 # ## DROT_TEST tests DROT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 06 September 2016 # # Author: # # John Burkardt # import numpy as np import platform n = 6 x = np.zeros ( n ) y = np.zeros ( n ) for i in range ( 0, n ): x[i] = float ( i + 1 ) y[i] = x[i] * x[i] - 12.0 print ( '' ) print ( 'DROT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' DROT carries out a Givens rotation.' ) print ( '' ) print ( ' Vectors X and Y' ) print ( '' ) for i in range ( 0, n ): print ( ' %6d %14.6g %14.6g' % ( i, x[i], y[i] ) ) c = 0.5 s = np.sqrt ( 1.0 - c * c ) xr, yr = drot ( n, x, 1, y, 1, c, s ) print ( '' ) print ( ' xr, yr = drot ( n, x, 1, y, 1, %g, %g )' % ( c, s ) ) print ( '' ) print ( ' Rotated vectors XR and YR' ) print ( '' ) for i in range ( 0, n ): print ( ' %6d %14.6g %14.6g' % ( i, xr[i], yr[i] ) ) t = np.arctan2 ( y[0], x[0] ) c = np.cos ( t ) s = np.sin ( t ) xr, yr = drot ( n, x, 1, y, 1, c, s ) print ( '' ) print ( ' xr, yr = drot ( n, x, 1, y, 1, %g, %g )' % ( c, s ) ) print ( '' ) print ( ' Rotated vectors XR and YR' ) print ( '' ) for i in range ( 0, n ): print ( ' %6d %14.6g %14.6g' % ( i, xr[i], yr[i] ) ) # # Terminate. # print ( '' ) print ( 'DROT_TEST' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) drot_test ( ) timestamp ( )