#! /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 ( )