1 SUBROUTINE ZLARTG( F, G, CS, SN, R ) 2* 3* -- LAPACK auxiliary routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* June 30, 1999 7* 8* .. Scalar Arguments .. 9 DOUBLE PRECISION CS 10 COMPLEX*16 F, G, R, SN 11* .. 12* 13* Purpose 14* ======= 15* 16* ZLARTG generates a plane rotation so that 17* 18* [ CS SN ] [ F ] [ R ] 19* [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1. 20* [ -SN CS ] [ G ] [ 0 ] 21* 22* This is a faster version of the BLAS1 routine ZROTG, except for 23* the following differences: 24* F and G are unchanged on return. 25* If G=0, then CS=1 and SN=0. 26* If F=0, then CS=0 and SN is chosen so that R is real. 27* 28* Arguments 29* ========= 30* 31* F (input) COMPLEX*16 32* The first component of vector to be rotated. 33* 34* G (input) COMPLEX*16 35* The second component of vector to be rotated. 36* 37* CS (output) DOUBLE PRECISION 38* The cosine of the rotation. 39* 40* SN (output) COMPLEX*16 41* The sine of the rotation. 42* 43* R (output) COMPLEX*16 44* The nonzero component of the rotated vector. 45* 46* Further Details 47* ======= ======= 48* 49* 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel 50* 51* ===================================================================== 52* 53* .. Parameters .. 54 DOUBLE PRECISION TWO, ONE, ZERO 55 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) 56 COMPLEX*16 CZERO 57 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 58* .. 59* .. Local Scalars .. 60 LOGICAL FIRST 61 INTEGER COUNT, I 62 DOUBLE PRECISION D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN, 63 $ SAFMN2, SAFMX2, SCALE 64 COMPLEX*16 FF, FS, GS 65* .. 66* .. External Functions .. 67 DOUBLE PRECISION DLAMCH, DLAPY2 68 EXTERNAL DLAMCH, DLAPY2 69* .. 70* .. Intrinsic Functions .. 71 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG, 72 $ MAX, SQRT 73* .. 74* .. Statement Functions .. 75 DOUBLE PRECISION ABS1, ABSSQ 76* .. 77* .. Save statement .. 78 SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 79* .. 80* .. Data statements .. 81 DATA FIRST / .TRUE. / 82* .. 83* .. Statement Function definitions .. 84 ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) ) 85 ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2 86* .. 87* .. Executable Statements .. 88* 89 IF( FIRST ) THEN 90 FIRST = .FALSE. 91 SAFMIN = DLAMCH( 'S' ) 92 EPS = DLAMCH( 'E' ) 93 SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / 94 $ LOG( DLAMCH( 'B' ) ) / TWO ) 95 SAFMX2 = ONE / SAFMN2 96 END IF 97 SCALE = MAX( ABS1( F ), ABS1( G ) ) 98 FS = F 99 GS = G 100 COUNT = 0 101 IF( SCALE.GE.SAFMX2 ) THEN 102 10 CONTINUE 103 COUNT = COUNT + 1 104 FS = FS*SAFMN2 105 GS = GS*SAFMN2 106 SCALE = SCALE*SAFMN2 107 IF( SCALE.GE.SAFMX2 ) 108 $ GO TO 10 109 ELSE IF( SCALE.LE.SAFMN2 ) THEN 110 IF( G.EQ.CZERO ) THEN 111 CS = ONE 112 SN = CZERO 113 R = F 114 RETURN 115 END IF 116 20 CONTINUE 117 COUNT = COUNT - 1 118 FS = FS*SAFMX2 119 GS = GS*SAFMX2 120 SCALE = SCALE*SAFMX2 121 IF( SCALE.LE.SAFMN2 ) 122 $ GO TO 20 123 END IF 124 F2 = ABSSQ( FS ) 125 G2 = ABSSQ( GS ) 126 IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN 127* 128* This is a rare case: F is very small. 129* 130 IF( F.EQ.CZERO ) THEN 131 CS = ZERO 132 R = DLAPY2( DBLE( G ), DIMAG( G ) ) 133* Do complex/real division explicitly with two real divisions 134 D = DLAPY2( DBLE( GS ), DIMAG( GS ) ) 135 SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D ) 136 RETURN 137 END IF 138 F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) ) 139* G2 and G2S are accurate 140* G2 is at least SAFMIN, and G2S is at least SAFMN2 141 G2S = SQRT( G2 ) 142* Error in CS from underflow in F2S is at most 143* UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS 144* If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN, 145* and so CS .lt. sqrt(SAFMIN) 146* If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN 147* and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS) 148* Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S 149 CS = F2S / G2S 150* Make sure abs(FF) = 1 151* Do complex/real division explicitly with 2 real divisions 152 IF( ABS1( F ).GT.ONE ) THEN 153 D = DLAPY2( DBLE( F ), DIMAG( F ) ) 154 FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D ) 155 ELSE 156 DR = SAFMX2*DBLE( F ) 157 DI = SAFMX2*DIMAG( F ) 158 D = DLAPY2( DR, DI ) 159 FF = DCMPLX( DR / D, DI / D ) 160 END IF 161 SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S ) 162 R = CS*F + SN*G 163 ELSE 164* 165* This is the most common case. 166* Neither F2 nor F2/G2 are less than SAFMIN 167* F2S cannot overflow, and it is accurate 168* 169 F2S = SQRT( ONE+G2 / F2 ) 170* Do the F2S(real)*FS(complex) multiply with two real multiplies 171 R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) ) 172 CS = ONE / F2S 173 D = F2 + G2 174* Do complex/real division explicitly with two real divisions 175 SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D ) 176 SN = SN*DCONJG( GS ) 177 IF( COUNT.NE.0 ) THEN 178 IF( COUNT.GT.0 ) THEN 179 DO 30 I = 1, COUNT 180 R = R*SAFMX2 181 30 CONTINUE 182 ELSE 183 DO 40 I = 1, -COUNT 184 R = R*SAFMN2 185 40 CONTINUE 186 END IF 187 END IF 188 END IF 189 RETURN 190* 191* End of ZLARTG 192* 193 END 194