[n,cx,incx,cy,incy,c,s]=csrot(n,cx,incx,cy,incy,c,s);

function [n,cx,incx,cy,incy,c,s]=csrot(n,cx,incx,cy,incy,c,s);
%***BEGIN PROLOGUE  CSROT
%***PURPOSE  Apply a plane Givens rotation.
%***LIBRARY   SLATEC (BLAS)
%***CATEGORY  D1B10
%***TYPE      COMPLEX (SROT-S, DROT-D, CSROT-C)
%***KEYWORDS  BLAS, GIVENS ROTATION, GIVENS TRANSFORMATION,
%             LINEAR ALGEBRA, PLANE ROTATION, VECTOR
%***AUTHOR  Dongarra, J., (ANL)
%***DESCRIPTION
%
%     CSROT applies the complex Givens rotation
%
%          (X)   ( C S,X)
%          (Y) = (-S C,Y)
%
%     N times where for I = 0,...,N-1
%
%          X = CX(LX+I*INCX)
%          Y = CY(LY+I*INCY),
%
%     where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
%     defined in a similar way using INCY.
%
%     Argument Description
%
%        N      (integer)  number of elements in each vector
%
%        CX     (complex array)  beginning of one vector
%
%        INCX   (integer)  memory spacing of successive elements
%               of vector CX
%
%        CY     (complex array)  beginning of the other vector
%
%        INCY   (integer)  memory spacing of successive elements
%               of vector CY
%
%        C      (real)  cosine term of the rotation
%
%        S      (real)  sine term of the rotation.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   810223  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920310  Corrected definition of LX in DESCRIPTION.  (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  CSROT
persistent ctemp i ix iy ; 

cx_shape=size(cx);cx=reshape(cx,1,[]);
cy_shape=size(cy);cy=reshape(cy,1,[]);
if isempty(ctemp), ctemp=0; end;
if isempty(i), i=0; end;
if isempty(ix), ix=0; end;
if isempty(iy), iy=0; end;
%***FIRST EXECUTABLE STATEMENT  CSROT
if( n<=0 )
cx_shape=zeros(cx_shape);cx_shape(:)=cx(1:numel(cx_shape));cx=cx_shape;
cy_shape=zeros(cy_shape);cy_shape(:)=cy(1:numel(cy_shape));cy=cy_shape;
return;
end;
if( incx==1 && incy==1 )
%
%     Code for both increments equal to 1.
%
for i = 1 : n;
ctemp = c.*cx(i) + s.*cy(i);
cy(i) = c.*cy(i) - s.*cx(i);
cx(i) = ctemp;
end; i = fix(n+1);
else;
%
%     Code for unequal increments or equal increments not equal to 1.
%
ix = 1;
iy = 1;
if( incx<0 )
ix =fix((-n+1).*incx + 1);
end;
if( incy<0 )
iy =fix((-n+1).*incy + 1);
end;
for i = 1 : n;
ctemp = c.*cx(ix) + s.*cy(iy);
cy(iy) = c.*cy(iy) - s.*cx(ix);
cx(ix) = ctemp;
ix = fix(ix + incx);
iy = fix(iy + incy);
end; i = fix(n+1);
cx_shape=zeros(cx_shape);cx_shape(:)=cx(1:numel(cx_shape));cx=cx_shape;
cy_shape=zeros(cy_shape);cy_shape(:)=cy(1:numel(cy_shape));cy=cy_shape;
return;
end;
cx_shape=zeros(cx_shape);cx_shape(:)=cx(1:numel(cx_shape));cx=cx_shape;
cy_shape=zeros(cy_shape);cy_shape(:)=cy(1:numel(cy_shape));cy=cy_shape;
end
%DECK CSSCAL