LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlargv.f
Go to the documentation of this file.
00001 *> \brief \b ZLARGV
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZLARGV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlargv.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlargv.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlargv.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZLARGV( N, X, INCX, Y, INCY, C, INCC )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INCC, INCX, INCY, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       DOUBLE PRECISION   C( * )
00028 *       COMPLEX*16         X( * ), Y( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> ZLARGV generates a vector of complex plane rotations with real
00038 *> cosines, determined by elements of the complex vectors x and y.
00039 *> For i = 1,2,...,n
00040 *>
00041 *>    (        c(i)   s(i) ) ( x(i) ) = ( r(i) )
00042 *>    ( -conjg(s(i))  c(i) ) ( y(i) ) = (   0  )
00043 *>
00044 *>    where c(i)**2 + ABS(s(i))**2 = 1
00045 *>
00046 *> The following conventions are used (these are the same as in ZLARTG,
00047 *> but differ from the BLAS1 routine ZROTG):
00048 *>    If y(i)=0, then c(i)=1 and s(i)=0.
00049 *>    If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.
00050 *> \endverbatim
00051 *
00052 *  Arguments:
00053 *  ==========
00054 *
00055 *> \param[in] N
00056 *> \verbatim
00057 *>          N is INTEGER
00058 *>          The number of plane rotations to be generated.
00059 *> \endverbatim
00060 *>
00061 *> \param[in,out] X
00062 *> \verbatim
00063 *>          X is COMPLEX*16 array, dimension (1+(N-1)*INCX)
00064 *>          On entry, the vector x.
00065 *>          On exit, x(i) is overwritten by r(i), for i = 1,...,n.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] INCX
00069 *> \verbatim
00070 *>          INCX is INTEGER
00071 *>          The increment between elements of X. INCX > 0.
00072 *> \endverbatim
00073 *>
00074 *> \param[in,out] Y
00075 *> \verbatim
00076 *>          Y is COMPLEX*16 array, dimension (1+(N-1)*INCY)
00077 *>          On entry, the vector y.
00078 *>          On exit, the sines of the plane rotations.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] INCY
00082 *> \verbatim
00083 *>          INCY is INTEGER
00084 *>          The increment between elements of Y. INCY > 0.
00085 *> \endverbatim
00086 *>
00087 *> \param[out] C
00088 *> \verbatim
00089 *>          C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
00090 *>          The cosines of the plane rotations.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] INCC
00094 *> \verbatim
00095 *>          INCC is INTEGER
00096 *>          The increment between elements of C. INCC > 0.
00097 *> \endverbatim
00098 *
00099 *  Authors:
00100 *  ========
00101 *
00102 *> \author Univ. of Tennessee 
00103 *> \author Univ. of California Berkeley 
00104 *> \author Univ. of Colorado Denver 
00105 *> \author NAG Ltd. 
00106 *
00107 *> \date November 2011
00108 *
00109 *> \ingroup complex16OTHERauxiliary
00110 *
00111 *> \par Further Details:
00112 *  =====================
00113 *>
00114 *> \verbatim
00115 *>
00116 *>  6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel
00117 *>
00118 *>  This version has a few statements commented out for thread safety
00119 *>  (machine parameters are computed on each entry). 10 feb 03, SJH.
00120 *> \endverbatim
00121 *>
00122 *  =====================================================================
00123       SUBROUTINE ZLARGV( N, X, INCX, Y, INCY, C, INCC )
00124 *
00125 *  -- LAPACK auxiliary routine (version 3.4.0) --
00126 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00127 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00128 *     November 2011
00129 *
00130 *     .. Scalar Arguments ..
00131       INTEGER            INCC, INCX, INCY, N
00132 *     ..
00133 *     .. Array Arguments ..
00134       DOUBLE PRECISION   C( * )
00135       COMPLEX*16         X( * ), Y( * )
00136 *     ..
00137 *
00138 *  =====================================================================
00139 *
00140 *     .. Parameters ..
00141       DOUBLE PRECISION   TWO, ONE, ZERO
00142       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
00143       COMPLEX*16         CZERO
00144       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
00145 *     ..
00146 *     .. Local Scalars ..
00147 *     LOGICAL            FIRST
00148 
00149       INTEGER            COUNT, I, IC, IX, IY, J
00150       DOUBLE PRECISION   CS, D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
00151      $                   SAFMN2, SAFMX2, SCALE
00152       COMPLEX*16         F, FF, FS, G, GS, R, SN
00153 *     ..
00154 *     .. External Functions ..
00155       DOUBLE PRECISION   DLAMCH, DLAPY2
00156       EXTERNAL           DLAMCH, DLAPY2
00157 *     ..
00158 *     .. Intrinsic Functions ..
00159       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
00160      $                   MAX, SQRT
00161 *     ..
00162 *     .. Statement Functions ..
00163       DOUBLE PRECISION   ABS1, ABSSQ
00164 *     ..
00165 *     .. Save statement ..
00166 *     SAVE               FIRST, SAFMX2, SAFMIN, SAFMN2
00167 *     ..
00168 *     .. Data statements ..
00169 *     DATA               FIRST / .TRUE. /
00170 *     ..
00171 *     .. Statement Function definitions ..
00172       ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
00173       ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
00174 *     ..
00175 *     .. Executable Statements ..
00176 *
00177 *     IF( FIRST ) THEN
00178 *        FIRST = .FALSE.
00179          SAFMIN = DLAMCH( 'S' )
00180          EPS = DLAMCH( 'E' )
00181          SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
00182      $            LOG( DLAMCH( 'B' ) ) / TWO )
00183          SAFMX2 = ONE / SAFMN2
00184 *     END IF
00185       IX = 1
00186       IY = 1
00187       IC = 1
00188       DO 60 I = 1, N
00189          F = X( IX )
00190          G = Y( IY )
00191 *
00192 *        Use identical algorithm as in ZLARTG
00193 *
00194          SCALE = MAX( ABS1( F ), ABS1( G ) )
00195          FS = F
00196          GS = G
00197          COUNT = 0
00198          IF( SCALE.GE.SAFMX2 ) THEN
00199    10       CONTINUE
00200             COUNT = COUNT + 1
00201             FS = FS*SAFMN2
00202             GS = GS*SAFMN2
00203             SCALE = SCALE*SAFMN2
00204             IF( SCALE.GE.SAFMX2 )
00205      $         GO TO 10
00206          ELSE IF( SCALE.LE.SAFMN2 ) THEN
00207             IF( G.EQ.CZERO ) THEN
00208                CS = ONE
00209                SN = CZERO
00210                R = F
00211                GO TO 50
00212             END IF
00213    20       CONTINUE
00214             COUNT = COUNT - 1
00215             FS = FS*SAFMX2
00216             GS = GS*SAFMX2
00217             SCALE = SCALE*SAFMX2
00218             IF( SCALE.LE.SAFMN2 )
00219      $         GO TO 20
00220          END IF
00221          F2 = ABSSQ( FS )
00222          G2 = ABSSQ( GS )
00223          IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
00224 *
00225 *           This is a rare case: F is very small.
00226 *
00227             IF( F.EQ.CZERO ) THEN
00228                CS = ZERO
00229                R = DLAPY2( DBLE( G ), DIMAG( G ) )
00230 *              Do complex/real division explicitly with two real
00231 *              divisions
00232                D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
00233                SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
00234                GO TO 50
00235             END IF
00236             F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
00237 *           G2 and G2S are accurate
00238 *           G2 is at least SAFMIN, and G2S is at least SAFMN2
00239             G2S = SQRT( G2 )
00240 *           Error in CS from underflow in F2S is at most
00241 *           UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
00242 *           If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
00243 *           and so CS .lt. sqrt(SAFMIN)
00244 *           If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
00245 *           and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
00246 *           Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
00247             CS = F2S / G2S
00248 *           Make sure abs(FF) = 1
00249 *           Do complex/real division explicitly with 2 real divisions
00250             IF( ABS1( F ).GT.ONE ) THEN
00251                D = DLAPY2( DBLE( F ), DIMAG( F ) )
00252                FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
00253             ELSE
00254                DR = SAFMX2*DBLE( F )
00255                DI = SAFMX2*DIMAG( F )
00256                D = DLAPY2( DR, DI )
00257                FF = DCMPLX( DR / D, DI / D )
00258             END IF
00259             SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
00260             R = CS*F + SN*G
00261          ELSE
00262 *
00263 *           This is the most common case.
00264 *           Neither F2 nor F2/G2 are less than SAFMIN
00265 *           F2S cannot overflow, and it is accurate
00266 *
00267             F2S = SQRT( ONE+G2 / F2 )
00268 *           Do the F2S(real)*FS(complex) multiply with two real
00269 *           multiplies
00270             R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
00271             CS = ONE / F2S
00272             D = F2 + G2
00273 *           Do complex/real division explicitly with two real divisions
00274             SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
00275             SN = SN*DCONJG( GS )
00276             IF( COUNT.NE.0 ) THEN
00277                IF( COUNT.GT.0 ) THEN
00278                   DO 30 J = 1, COUNT
00279                      R = R*SAFMX2
00280    30             CONTINUE
00281                ELSE
00282                   DO 40 J = 1, -COUNT
00283                      R = R*SAFMN2
00284    40             CONTINUE
00285                END IF
00286             END IF
00287          END IF
00288    50    CONTINUE
00289          C( IC ) = CS
00290          Y( IY ) = SN
00291          X( IX ) = R
00292          IC = IC + INCC
00293          IY = IY + INCY
00294          IX = IX + INCX
00295    60 CONTINUE
00296       RETURN
00297 *
00298 *     End of ZLARGV
00299 *
00300       END
 All Files Functions