de/df7/zlargv_8f_source.html

*> \brief \b ZLARGV generates a vector of plane rotations with real cosines and complex sines.

*

*  =========== DOCUMENTATION ===========

*

* Online html documentation available at

*            http://www.netlib.org/lapack/explore-html/

*

*> \htmlonly

*> Download ZLARGV + dependencies

*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlargv.f">

*> [TGZ]</a>

*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlargv.f">

*> [ZIP]</a>

*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlargv.f">

*> [TXT]</a>

*> \endhtmlonly

*

*  Definition:

*  ===========

*

*       SUBROUTINE ZLARGV( N, X, INCX, Y, INCY, C, INCC )

*

*       .. Scalar Arguments ..

*       INTEGER            INCC, INCX, INCY, N

*       ..

*       .. Array Arguments ..

*       DOUBLE PRECISION   C( * )

*       COMPLEX*16         X( * ), Y( * )

*       ..

*

*

*> \par Purpose:

*  =============

*>

*> \verbatim

*>

*> ZLARGV generates a vector of complex plane rotations with real

*> cosines, determined by elements of the complex vectors x and y.

*> For i = 1,2,...,n

*>

*>    (        c(i)   s(i) ) ( x(i) ) = ( r(i) )

*>    ( -conjg(s(i))  c(i) ) ( y(i) ) = (   0  )

*>

*>    where c(i)**2 + ABS(s(i))**2 = 1

*>

*> The following conventions are used (these are the same as in ZLARTG,

*> but differ from the BLAS1 routine ZROTG):

*>    If y(i)=0, then c(i)=1 and s(i)=0.

*>    If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.

*> \endverbatim

*

*  Arguments:

*  ==========

*

*> \param[in] N

*> \verbatim

*>          N is INTEGER

*>          The number of plane rotations to be generated.

*> \endverbatim

*>

*> \param[in,out] X

*> \verbatim

*>          X is COMPLEX*16 array, dimension (1+(N-1)*INCX)

*>          On entry, the vector x.

*>          On exit, x(i) is overwritten by r(i), for i = 1,...,n.

*> \endverbatim

*>

*> \param[in] INCX

*> \verbatim

*>          INCX is INTEGER

*>          The increment between elements of X. INCX > 0.

*> \endverbatim

*>

*> \param[in,out] Y

*> \verbatim

*>          Y is COMPLEX*16 array, dimension (1+(N-1)*INCY)

*>          On entry, the vector y.

*>          On exit, the sines of the plane rotations.

*> \endverbatim

*>

*> \param[in] INCY

*> \verbatim

*>          INCY is INTEGER

*>          The increment between elements of Y. INCY > 0.

*> \endverbatim

*>

*> \param[out] C

*> \verbatim

*>          C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)

*>          The cosines of the plane rotations.

*> \endverbatim

*>

*> \param[in] INCC

*> \verbatim

*>          INCC is INTEGER

*>          The increment between elements of C. INCC > 0.

*> \endverbatim

*

*  Authors:

*  ========

*

*> \author Univ. of Tennessee

*> \author Univ. of California Berkeley

*> \author Univ. of Colorado Denver

*> \author NAG Ltd.

*

*> \ingroup largv

*

*> \par Further Details:

*  =====================

*>

*> \verbatim

*>

*>  6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel

*>

*>  This version has a few statements commented out for thread safety

*>  (machine parameters are computed on each entry). 10 feb 03, SJH.

*> \endverbatim

*>

*  =====================================================================

      SUBROUTINE zlargv( N, X, INCX, Y, INCY, C, INCC )

*

*  -- LAPACK auxiliary routine --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*

*     .. Scalar Arguments ..

      INTEGER            INCC, INCX, INCY, N

*     ..

*     .. Array Arguments ..

      DOUBLE PRECISION   C( * )

      COMPLEX*16         X( * ), Y( * )

*     ..

*

*  =====================================================================

*

*     .. Parameters ..

      DOUBLE PRECISION   TWO, ONE, ZERO

      parameter( two = 2.0d+0, one = 1.0d+0, zero = 0.0d+0 )

      COMPLEX*16         CZERO

      parameter( czero = ( 0.0d+0, 0.0d+0 ) )

*     ..

*     .. Local Scalars ..

*     LOGICAL            FIRST


      INTEGER            COUNT, I, IC, IX, IY, J

      DOUBLE PRECISION   CS, D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,

     $                   SAFMN2, SAFMX2, SCALE

      COMPLEX*16         F, FF, FS, G, GS, R, SN

*     ..

*     .. External Functions ..

      DOUBLE PRECISION   DLAMCH, DLAPY2

      EXTERNAL           dlamch, dlapy2

*     ..

*     .. Intrinsic Functions ..

      INTRINSIC          abs, dble, dcmplx, dconjg, dimag, int, log,

     $                   max, sqrt

*     ..

*     .. Statement Functions ..

      DOUBLE PRECISION   ABS1, ABSSQ

*     ..

*     .. Save statement ..

*     SAVE               FIRST, SAFMX2, SAFMIN, SAFMN2

*     ..

*     .. Data statements ..

*     DATA               FIRST / .TRUE. /

*     ..

*     .. Statement Function definitions ..

      abs1( ff ) = max( abs( dble( ff ) ), abs( dimag( ff ) ) )

      abssq( ff ) = dble( ff )**2 + dimag( ff )**2

*     ..

*     .. Executable Statements ..

*

*     IF( FIRST ) THEN

*        FIRST = .FALSE.

         safmin = dlamch( 'S' )

         eps = dlamch( 'E' )

         safmn2 = dlamch( 'B' )**int( log( safmin / eps ) /

     $            log( dlamch( 'B' ) ) / two )

         safmx2 = one / safmn2

*     END IF

      ix = 1

      iy = 1

      ic = 1

      DO 60 i = 1, n

         f = x( ix )

         g = y( iy )

*

*        Use identical algorithm as in ZLARTG

*

         scale = max( abs1( f ), abs1( g ) )

         fs = f

         gs = g

         count = 0

         IF( scale.GE.safmx2 ) THEN

   10       CONTINUE

            count = count + 1

            fs = fs*safmn2

            gs = gs*safmn2

            scale = scale*safmn2

            IF( scale.GE.safmx2 .AND. count .LT. 20 )

     $         GO TO 10

         ELSE IF( scale.LE.safmn2 ) THEN

            IF( g.EQ.czero ) THEN

               cs = one

               sn = czero

               r = f

               GO TO 50

            END IF

   20       CONTINUE

            count = count - 1

            fs = fs*safmx2

            gs = gs*safmx2

            scale = scale*safmx2

            IF( scale.LE.safmn2 )

     $         GO TO 20

         END IF

         f2 = abssq( fs )

         g2 = abssq( gs )

         IF( f2.LE.max( g2, one )*safmin ) THEN

*

*           This is a rare case: F is very small.

*

            IF( f.EQ.czero ) THEN

               cs = zero

               r = dlapy2( dble( g ), dimag( g ) )

*              Do complex/real division explicitly with two real

*              divisions

               d = dlapy2( dble( gs ), dimag( gs ) )

               sn = dcmplx( dble( gs ) / d, -dimag( gs ) / d )

               GO TO 50

            END IF

            f2s = dlapy2( dble( fs ), dimag( fs ) )

*           G2 and G2S are accurate

*           G2 is at least SAFMIN, and G2S is at least SAFMN2

            g2s = sqrt( g2 )

*           Error in CS from underflow in F2S is at most

*           UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS

*           If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,

*           and so CS .lt. sqrt(SAFMIN)

*           If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN

*           and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)

*           Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S

            cs = f2s / g2s

*           Make sure abs(FF) = 1

*           Do complex/real division explicitly with 2 real divisions

            IF( abs1( f ).GT.one ) THEN

               d = dlapy2( dble( f ), dimag( f ) )

               ff = dcmplx( dble( f ) / d, dimag( f ) / d )

            ELSE

               dr = safmx2*dble( f )

               di = safmx2*dimag( f )

               d = dlapy2( dr, di )

               ff = dcmplx( dr / d, di / d )

            END IF

            sn = ff*dcmplx( dble( gs ) / g2s, -dimag( gs ) / g2s )

            r = cs*f + sn*g

         ELSE

*

*           This is the most common case.

*           Neither F2 nor F2/G2 are less than SAFMIN

*           F2S cannot overflow, and it is accurate

*

            f2s = sqrt( one+g2 / f2 )

*           Do the F2S(real)*FS(complex) multiply with two real

*           multiplies

            r = dcmplx( f2s*dble( fs ), f2s*dimag( fs ) )

            cs = one / f2s

            d = f2 + g2

*           Do complex/real division explicitly with two real divisions

            sn = dcmplx( dble( r ) / d, dimag( r ) / d )

            sn = sn*dconjg( gs )

            IF( count.NE.0 ) THEN

               IF( count.GT.0 ) THEN

                  DO 30 j = 1, count

                     r = r*safmx2

   30             CONTINUE

               ELSE

                  DO 40 j = 1, -count

                     r = r*safmn2

   40             CONTINUE

               END IF

            END IF

         END IF

   50    CONTINUE

         c( ic ) = cs

         y( iy ) = sn

         x( ix ) = r

         ic = ic + incc

         iy = iy + incy

         ix = ix + incx

   60 CONTINUE

      RETURN

*

*     End of ZLARGV

*

      END

zlargv
subroutine zlargv(n, x, incx, y, incy, c, incc)
ZLARGV generates a vector of plane rotations with real cosines and complex sines.
Definition zlargv.f:122