LAPACK: SRC/zlargv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLARGV( N, X, INCX, Y, INCY, C, INCC )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INCC, INCX, INCY, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   C( * )
00013       COMPLEX*16         X( * ), Y( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZLARGV generates a vector of complex plane rotations with real
00020 *  cosines, determined by elements of the complex vectors x and y.
00021 *  For i = 1,2,...,n
00022 *
00023 *     (        c(i)   s(i) ) ( x(i) ) = ( r(i) )
00024 *     ( -conjg(s(i))  c(i) ) ( y(i) ) = (   0  )
00025 *
00026 *     where c(i)**2 + ABS(s(i))**2 = 1
00027 *
00028 *  The following conventions are used (these are the same as in ZLARTG,
00029 *  but differ from the BLAS1 routine ZROTG):
00030 *     If y(i)=0, then c(i)=1 and s(i)=0.
00031 *     If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.
00032 *
00033 *  Arguments
00034 *  =========
00035 *
00036 *  N       (input) INTEGER
00037 *          The number of plane rotations to be generated.
00038 *
00039 *  X       (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCX)
00040 *          On entry, the vector x.
00041 *          On exit, x(i) is overwritten by r(i), for i = 1,...,n.
00042 *
00043 *  INCX    (input) INTEGER
00044 *          The increment between elements of X. INCX > 0.
00045 *
00046 *  Y       (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCY)
00047 *          On entry, the vector y.
00048 *          On exit, the sines of the plane rotations.
00049 *
00050 *  INCY    (input) INTEGER
00051 *          The increment between elements of Y. INCY > 0.
00052 *
00053 *  C       (output) DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
00054 *          The cosines of the plane rotations.
00055 *
00056 *  INCC    (input) INTEGER
00057 *          The increment between elements of C. INCC > 0.
00058 *
00059 *  Further Details
00060 *  ======= =======
00061 *
00062 *  6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel
00063 *
00064 *  This version has a few statements commented out for thread safety
00065 *  (machine parameters are computed on each entry). 10 feb 03, SJH.
00066 *
00067 *  =====================================================================
00068 *
00069 *     .. Parameters ..
00070       DOUBLE PRECISION   TWO, ONE, ZERO
00071       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
00072       COMPLEX*16         CZERO
00073       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
00074 *     ..
00075 *     .. Local Scalars ..
00076 *     LOGICAL            FIRST
00077 
00078       INTEGER            COUNT, I, IC, IX, IY, J
00079       DOUBLE PRECISION   CS, D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
00080      $                   SAFMN2, SAFMX2, SCALE
00081       COMPLEX*16         F, FF, FS, G, GS, R, SN
00082 *     ..
00083 *     .. External Functions ..
00084       DOUBLE PRECISION   DLAMCH, DLAPY2
00085       EXTERNAL           DLAMCH, DLAPY2
00086 *     ..
00087 *     .. Intrinsic Functions ..
00088       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
00089      $                   MAX, SQRT
00090 *     ..
00091 *     .. Statement Functions ..
00092       DOUBLE PRECISION   ABS1, ABSSQ
00093 *     ..
00094 *     .. Save statement ..
00095 *     SAVE               FIRST, SAFMX2, SAFMIN, SAFMN2
00096 *     ..
00097 *     .. Data statements ..
00098 *     DATA               FIRST / .TRUE. /
00099 *     ..
00100 *     .. Statement Function definitions ..
00101       ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
00102       ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
00103 *     ..
00104 *     .. Executable Statements ..
00105 *
00106 *     IF( FIRST ) THEN
00107 *        FIRST = .FALSE.
00108          SAFMIN = DLAMCH( 'S' )
00109          EPS = DLAMCH( 'E' )
00110          SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
00111      $            LOG( DLAMCH( 'B' ) ) / TWO )
00112          SAFMX2 = ONE / SAFMN2
00113 *     END IF
00114       IX = 1
00115       IY = 1
00116       IC = 1
00117       DO 60 I = 1, N
00118          F = X( IX )
00119          G = Y( IY )
00120 *
00121 *        Use identical algorithm as in ZLARTG
00122 *
00123          SCALE = MAX( ABS1( F ), ABS1( G ) )
00124          FS = F
00125          GS = G
00126          COUNT = 0
00127          IF( SCALE.GE.SAFMX2 ) THEN
00128    10       CONTINUE
00129             COUNT = COUNT + 1
00130             FS = FS*SAFMN2
00131             GS = GS*SAFMN2
00132             SCALE = SCALE*SAFMN2
00133             IF( SCALE.GE.SAFMX2 )
00134      $         GO TO 10
00135          ELSE IF( SCALE.LE.SAFMN2 ) THEN
00136             IF( G.EQ.CZERO ) THEN
00137                CS = ONE
00138                SN = CZERO
00139                R = F
00140                GO TO 50
00141             END IF
00142    20       CONTINUE
00143             COUNT = COUNT - 1
00144             FS = FS*SAFMX2
00145             GS = GS*SAFMX2
00146             SCALE = SCALE*SAFMX2
00147             IF( SCALE.LE.SAFMN2 )
00148      $         GO TO 20
00149          END IF
00150          F2 = ABSSQ( FS )
00151          G2 = ABSSQ( GS )
00152          IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
00153 *
00154 *           This is a rare case: F is very small.
00155 *
00156             IF( F.EQ.CZERO ) THEN
00157                CS = ZERO
00158                R = DLAPY2( DBLE( G ), DIMAG( G ) )
00159 *              Do complex/real division explicitly with two real
00160 *              divisions
00161                D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
00162                SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
00163                GO TO 50
00164             END IF
00165             F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
00166 *           G2 and G2S are accurate
00167 *           G2 is at least SAFMIN, and G2S is at least SAFMN2
00168             G2S = SQRT( G2 )
00169 *           Error in CS from underflow in F2S is at most
00170 *           UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
00171 *           If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
00172 *           and so CS .lt. sqrt(SAFMIN)
00173 *           If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
00174 *           and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
00175 *           Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
00176             CS = F2S / G2S
00177 *           Make sure abs(FF) = 1
00178 *           Do complex/real division explicitly with 2 real divisions
00179             IF( ABS1( F ).GT.ONE ) THEN
00180                D = DLAPY2( DBLE( F ), DIMAG( F ) )
00181                FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
00182             ELSE
00183                DR = SAFMX2*DBLE( F )
00184                DI = SAFMX2*DIMAG( F )
00185                D = DLAPY2( DR, DI )
00186                FF = DCMPLX( DR / D, DI / D )
00187             END IF
00188             SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
00189             R = CS*F + SN*G
00190          ELSE
00191 *
00192 *           This is the most common case.
00193 *           Neither F2 nor F2/G2 are less than SAFMIN
00194 *           F2S cannot overflow, and it is accurate
00195 *
00196             F2S = SQRT( ONE+G2 / F2 )
00197 *           Do the F2S(real)*FS(complex) multiply with two real
00198 *           multiplies
00199             R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
00200             CS = ONE / F2S
00201             D = F2 + G2
00202 *           Do complex/real division explicitly with two real divisions
00203             SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
00204             SN = SN*DCONJG( GS )
00205             IF( COUNT.NE.0 ) THEN
00206                IF( COUNT.GT.0 ) THEN
00207                   DO 30 J = 1, COUNT
00208                      R = R*SAFMX2
00209    30             CONTINUE
00210                ELSE
00211                   DO 40 J = 1, -COUNT
00212                      R = R*SAFMN2
00213    40             CONTINUE
00214                END IF
00215             END IF
00216          END IF
00217    50    CONTINUE
00218          C( IC ) = CS
00219          Y( IY ) = SN
00220          X( IX ) = R
00221          IC = IC + INCC
00222          IY = IY + INCY
00223          IX = IX + INCX
00224    60 CONTINUE
00225       RETURN
00226 *
00227 *     End of ZLARGV
00228 *
00229       END