![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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