C     THIS FILE CONTAINS THE ROUTINES THAT CONSTRUCT AND 
C     APPLY GIVENS ROTATIONS IN THE FOUR ARITHMETICS
C     SROT, SROTG, DROT, DROTG ARE EXTRACTED FROM THE BLAS LIBRARY 
C     CROT, ZROT ARE EXTRACTED FROM THE LAPACK LIBRARY
C
      subroutine srot (n,sx,incx,sy,incy,c,s)
c
c     applies a plane rotation.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sx(*),sy(*),stemp,c,s
      integer i,incx,incy,ix,iy,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c     code for unequal increments or equal increments not equal
c     to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        stemp  = c*sx(ix) + s*sy(iy)
        sy(iy) = c*sy(iy) - s*sx(ix)
        sx(ix) = stemp
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c     code for both increments equal to 1
c
   20 do 30 i = 1,n
        stemp = c*sx(i) + s*sy(i)
        sy(i) = c*sy(i) - s*sx(i)
        sx(i) = stemp
   30 continue
      return
      end
C
C
      subroutine srotg(sa,sb,c,s)
c
c     construct givens plane rotation.
c     jack dongarra, linpack, 3/11/78.
c
      real sa,sb,c,s,roe,scale,r,z
c
      roe = sb
      if( abs(sa) .gt. abs(sb) ) roe = sa
      scale = abs(sa) + abs(sb)
C
      if( scale .ne. 0.0 ) go to 10
C
      c = 1.0
      s = 0.0
      r = 0.0
      z = 0.0
      go to 20
C
   10 r = scale*sqrt((sa/scale)**2 + (sb/scale)**2)
      r = sign(1.0,roe)*r
      c = sa/r
      s = sb/r
      z = 1.0
      if( abs(sa) .gt. abs(sb) ) z = s
      if( abs(sb) .ge. abs(sa) .and. c .ne. 0.0 ) z = 1.0/c
C
   20 sa = r
      sb = z
      return
      end
C
C
      subroutine  drot (n,dx,incx,dy,incy,c,s)
c
c     applies a plane rotation.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      double precision dx(*),dy(*),dtemp,c,s
      integer i,incx,incy,ix,iy,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c     code for unequal increments or equal increments not equal
c     to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = c*dx(ix) + s*dy(iy)
        dy(iy) = c*dy(iy) - s*dx(ix)
        dx(ix) = dtemp
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c     code for both increments equal to 1
c
   20 do 30 i = 1,n
        dtemp = c*dx(i) + s*dy(i)
        dy(i) = c*dy(i) - s*dx(i)
        dx(i) = dtemp
   30 continue
      return
      end
C
C
      subroutine drotg(da,db,c,s)
c
c     construct givens plane rotation.
c     jack dongarra, linpack, 3/11/78.
c
      double precision da,db,c,s,roe,scale,r,z
c
      roe = db
      if( dabs(da) .gt. dabs(db) ) roe = da
      scale = dabs(da) + dabs(db)
C
      if( scale .ne. 0.0d0 ) go to 10
C
      c = 1.0d0
      s = 0.0d0
      r = 0.0d0
      z = 0.0d0
      go to 20
C
   10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2)
      r = dsign(1.0d0,roe)*r
      c = da/r
      s = db/r
      z = 1.0d0
      if( dabs(da) .gt. dabs(db) ) z = s
      if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c
C
   20 da = r
      db = z
      return
      end
C
C
      subroutine crot( n, cx, incx, cy, incy, c, s )
C
C  -- lapack auxiliary routine (version 3.0) --
C     univ. of tennessee, univ. of california berkeley, nag ltd.,
C     courant institute, argonne national lab, and rice university
C     october 31, 1992
C
C     .. scalar arguments ..
      integer            incx, incy, n
      real               c
      complex            s
C     ..
C     .. array arguments ..
      complex            cx( * ), cy( * )
C     ..
C
C  purpose
C  =======
C
C  crot   applies a plane rotation, where the cos (c) is real and the
C  sin (s) is complex, and the vectors cx and cy are complex.
C
C  arguments
C  =========
C
C  n       (input) integer
C          the number of elements in the vectors cx and cy.
C
C  cx      (input/output) complex array, dimension (n)
C          on input, the vector x.
C          on output, cx is overwritten with c*x + s*y.
C
C  incx    (input) integer
C          the increment between successive values of cy.  incx <> 0.
C
C  cy      (input/output) complex array, dimension (n)
C          on input, the vector y.
C          on output, cy is overwritten with -conjg(s)*x + c*y.
C
C  incy    (input) integer
C          the increment between successive values of cy.  incx <> 0.
C
C  c       (input) real
C  s       (input) complex
C          c and s define a rotation
C             [  c          s  ]
C             [ -conjg(s)   c  ]
C          where c*c + s*conjg(s) = 1.0.
C
C =====================================================================
C
C     .. local scalars ..
      integer            i, ix, iy
      complex            stemp
C     ..
C     .. intrinsic functions ..
      intrinsic          conjg
C     ..
C     .. executable statements ..
C
      if( n.le.0 )return
      if( incx.eq.1 .and. incy.eq.1 )go to 20
C
C     code for unequal increments or equal increments not equal 
C     to 1
C
      ix = 1
      iy = 1
      if( incx.lt.0 )ix = ( -n+1 )*incx + 1
      if( incy.lt.0 )iy = ( -n+1 )*incy + 1
      do 10 i = 1, n
         stemp  = c*cx(ix) + s*cy(iy)
         cy(iy) = c*cy(iy) - conjg(s)*cx(ix)
         cx(ix) = stemp
         ix = ix + incx
         iy = iy + incy
   10 continue
      return
C
C     code for both increments equal to 1
C
   20 do 30 i = 1, n
         stemp = c*cx(i) + s*cy(i)
         cy(i) = c*cy(i) - conjg(s)*cx(i)
         cx(i) = stemp
   30 continue
      return
      end
C
C
      subroutine crotg(ca,cb,c,s)
      complex ca,cb,s
      real c
      real norm,scale
      complex alpha
C
      if (cabs(ca) .ne. 0.) go to 10
C
      c = 0.
      s = (1.,0.)
      ca = cb
      go to 20
C
   10 continue
      scale = cabs(ca) + cabs(cb)
      norm = scale * sqrt((cabs(ca/scale))**2 + (cabs(cb/scale))**2)
      alpha = ca /cabs(ca)
      c = cabs(ca) / norm
      s = alpha * conjg(cb) / norm
      ca = alpha * norm
C
   20 continue
      return
      end
C
C
      subroutine zrot( n, cx, incx, cy, incy, c, s )
C
C  -- lapack auxiliary routine (version 3.0) --
C     univ. of tennessee, univ. of california berkeley, nag ltd.,
C     courant institute, argonne national lab, and rice university
C     october 31, 1992
C
C     .. scalar arguments ..
      integer            incx, incy, n
      double precision   c
      complex*16         s
C     ..
C     .. array arguments ..
      complex*16         cx( * ), cy( * )
C     ..
C
C  purpose
C  =======
C
C  zrot   applies a plane rotation, where the cos (c) is real and the
C  sin (s) is complex, and the vectors cx and cy are complex.
C
C  arguments
C  =========
C
C  n       (input) integer
C          the number of elements in the vectors cx and cy.
C
C  cx      (input/output) complex*16 array, dimension (n)
C          on input, the vector x.
C          on output, cx is overwritten with c*x + s*y.
C
C  incx    (input) integer
C          the increment between successive values of cy.  incx <> 0.
C
C  cy      (input/output) complex*16 array, dimension (n)
C          on input, the vector y.
C          on output, cy is overwritten with -conjg(s)*x + c*y.
C
C  incy    (input) integer
C          the increment between successive values of cy.  incx <> 0.
C
C  c       (input) double precision
C  s       (input) complex*16
C          c and s define a rotation
C             [  c          s  ]
C             [ -conjg(s)   c  ]
C          where c*c + s*conjg(s) = 1.0.
C
C =====================================================================
C
C     .. local scalars ..
      integer            i, ix, iy
      complex*16         stemp
C     ..
C     .. intrinsic functions ..
      intrinsic          dconjg
C     ..
C     .. executable statements ..
C
      if( n.le.0 )return
      if( incx.eq.1 .and. incy.eq.1 )go to 20
C
C     code for unequal increments or equal increments not equal to 1
C
      ix = 1
      iy = 1
      if( incx.lt.0 )ix = ( -n+1 )*incx + 1
      if( incy.lt.0 )iy = ( -n+1 )*incy + 1
      do 10 i = 1, n
         stemp  = c*cx(ix) + s*cy(iy)
         cy(iy) = c*cy(iy) - dconjg(s)*cx(ix)
         cx(ix) = stemp
         ix = ix + incx
         iy = iy + incy
   10 continue
      return
C
C     code for both increments equal to 1
C
   20 continue
      do 30 i = 1, n
         stemp = c*cx(i) + s*cy(i)
         cy(i) = c*cy(i) - dconjg(s)*cx(i)
         cx(i) = stemp
   30 continue
      return
      end
C
C
      subroutine zrotg(ca,cb,c,s)
      double complex ca,cb,s
      double precision c
      double precision norm,scale
      double complex alpha
C
      if (cdabs(ca) .ne. 0.0d0) go to 10
C
      c = 0.0d0
      s = (1.0d0,0.0d0)
      ca = cb
      go to 20
C
   10 continue
      scale = cdabs(ca) + cdabs(cb)
      norm = scale*dsqrt((cdabs(ca/dcmplx(scale,0.0d0)))**2 +
     &                      (cdabs(cb/dcmplx(scale,0.0d0)))**2)
      alpha = ca /cdabs(ca)
      c = cdabs(ca) / norm
      s = alpha * dconjg(cb) / norm
      ca = alpha * norm
C
   20 continue
      return
      end