|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CLARGV( 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 REAL C( * ) 00013 COMPLEX X( * ), Y( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CLARGV 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 CLARTG, 00029 * but differ from the BLAS1 routine CROTG): 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 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 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) REAL 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 REAL TWO, ONE, ZERO 00071 PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) 00072 COMPLEX CZERO 00073 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00074 * .. 00075 * .. Local Scalars .. 00076 * LOGICAL FIRST 00077 INTEGER COUNT, I, IC, IX, IY, J 00078 REAL CS, D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN, 00079 $ SAFMN2, SAFMX2, SCALE 00080 COMPLEX F, FF, FS, G, GS, R, SN 00081 * .. 00082 * .. External Functions .. 00083 REAL SLAMCH, SLAPY2 00084 EXTERNAL SLAMCH, SLAPY2 00085 * .. 00086 * .. Intrinsic Functions .. 00087 INTRINSIC ABS, AIMAG, CMPLX, CONJG, INT, LOG, MAX, REAL, 00088 $ SQRT 00089 * .. 00090 * .. Statement Functions .. 00091 REAL ABS1, ABSSQ 00092 * .. 00093 * .. Save statement .. 00094 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 00095 * .. 00096 * .. Data statements .. 00097 * DATA FIRST / .TRUE. / 00098 * .. 00099 * .. Statement Function definitions .. 00100 ABS1( FF ) = MAX( ABS( REAL( FF ) ), ABS( AIMAG( FF ) ) ) 00101 ABSSQ( FF ) = REAL( FF )**2 + AIMAG( FF )**2 00102 * .. 00103 * .. Executable Statements .. 00104 * 00105 * IF( FIRST ) THEN 00106 * FIRST = .FALSE. 00107 SAFMIN = SLAMCH( 'S' ) 00108 EPS = SLAMCH( 'E' ) 00109 SAFMN2 = SLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / 00110 $ LOG( SLAMCH( 'B' ) ) / TWO ) 00111 SAFMX2 = ONE / SAFMN2 00112 * END IF 00113 IX = 1 00114 IY = 1 00115 IC = 1 00116 DO 60 I = 1, N 00117 F = X( IX ) 00118 G = Y( IY ) 00119 * 00120 * Use identical algorithm as in CLARTG 00121 * 00122 SCALE = MAX( ABS1( F ), ABS1( G ) ) 00123 FS = F 00124 GS = G 00125 COUNT = 0 00126 IF( SCALE.GE.SAFMX2 ) THEN 00127 10 CONTINUE 00128 COUNT = COUNT + 1 00129 FS = FS*SAFMN2 00130 GS = GS*SAFMN2 00131 SCALE = SCALE*SAFMN2 00132 IF( SCALE.GE.SAFMX2 ) 00133 $ GO TO 10 00134 ELSE IF( SCALE.LE.SAFMN2 ) THEN 00135 IF( G.EQ.CZERO ) THEN 00136 CS = ONE 00137 SN = CZERO 00138 R = F 00139 GO TO 50 00140 END IF 00141 20 CONTINUE 00142 COUNT = COUNT - 1 00143 FS = FS*SAFMX2 00144 GS = GS*SAFMX2 00145 SCALE = SCALE*SAFMX2 00146 IF( SCALE.LE.SAFMN2 ) 00147 $ GO TO 20 00148 END IF 00149 F2 = ABSSQ( FS ) 00150 G2 = ABSSQ( GS ) 00151 IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN 00152 * 00153 * This is a rare case: F is very small. 00154 * 00155 IF( F.EQ.CZERO ) THEN 00156 CS = ZERO 00157 R = SLAPY2( REAL( G ), AIMAG( G ) ) 00158 * Do complex/real division explicitly with two real 00159 * divisions 00160 D = SLAPY2( REAL( GS ), AIMAG( GS ) ) 00161 SN = CMPLX( REAL( GS ) / D, -AIMAG( GS ) / D ) 00162 GO TO 50 00163 END IF 00164 F2S = SLAPY2( REAL( FS ), AIMAG( FS ) ) 00165 * G2 and G2S are accurate 00166 * G2 is at least SAFMIN, and G2S is at least SAFMN2 00167 G2S = SQRT( G2 ) 00168 * Error in CS from underflow in F2S is at most 00169 * UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS 00170 * If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN, 00171 * and so CS .lt. sqrt(SAFMIN) 00172 * If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN 00173 * and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS) 00174 * Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S 00175 CS = F2S / G2S 00176 * Make sure abs(FF) = 1 00177 * Do complex/real division explicitly with 2 real divisions 00178 IF( ABS1( F ).GT.ONE ) THEN 00179 D = SLAPY2( REAL( F ), AIMAG( F ) ) 00180 FF = CMPLX( REAL( F ) / D, AIMAG( F ) / D ) 00181 ELSE 00182 DR = SAFMX2*REAL( F ) 00183 DI = SAFMX2*AIMAG( F ) 00184 D = SLAPY2( DR, DI ) 00185 FF = CMPLX( DR / D, DI / D ) 00186 END IF 00187 SN = FF*CMPLX( REAL( GS ) / G2S, -AIMAG( GS ) / G2S ) 00188 R = CS*F + SN*G 00189 ELSE 00190 * 00191 * This is the most common case. 00192 * Neither F2 nor F2/G2 are less than SAFMIN 00193 * F2S cannot overflow, and it is accurate 00194 * 00195 F2S = SQRT( ONE+G2 / F2 ) 00196 * Do the F2S(real)*FS(complex) multiply with two real 00197 * multiplies 00198 R = CMPLX( F2S*REAL( FS ), F2S*AIMAG( FS ) ) 00199 CS = ONE / F2S 00200 D = F2 + G2 00201 * Do complex/real division explicitly with two real divisions 00202 SN = CMPLX( REAL( R ) / D, AIMAG( R ) / D ) 00203 SN = SN*CONJG( GS ) 00204 IF( COUNT.NE.0 ) THEN 00205 IF( COUNT.GT.0 ) THEN 00206 DO 30 J = 1, COUNT 00207 R = R*SAFMX2 00208 30 CONTINUE 00209 ELSE 00210 DO 40 J = 1, -COUNT 00211 R = R*SAFMN2 00212 40 CONTINUE 00213 END IF 00214 END IF 00215 END IF 00216 50 CONTINUE 00217 C( IC ) = CS 00218 Y( IY ) = SN 00219 X( IX ) = R 00220 IC = IC + INCC 00221 IY = IY + INCY 00222 IX = IX + INCX 00223 60 CONTINUE 00224 RETURN 00225 * 00226 * End of CLARGV 00227 * 00228 END
1.7.3 
