diff options
Diffstat (limited to 'src/fortran/lapack/zlartg.f')
| -rw-r--r-- | src/fortran/lapack/zlartg.f | 195 |
1 files changed, 195 insertions, 0 deletions
diff --git a/src/fortran/lapack/zlartg.f b/src/fortran/lapack/zlartg.f new file mode 100644 index 0000000..6d3a850 --- /dev/null +++ b/src/fortran/lapack/zlartg.f @@ -0,0 +1,195 @@ + SUBROUTINE ZLARTG( F, G, CS, SN, R ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + DOUBLE PRECISION CS + COMPLEX*16 F, G, R, SN +* .. +* +* Purpose +* ======= +* +* ZLARTG generates a plane rotation so that +* +* [ CS SN ] [ F ] [ R ] +* [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1. +* [ -SN CS ] [ G ] [ 0 ] +* +* This is a faster version of the BLAS1 routine ZROTG, except for +* the following differences: +* F and G are unchanged on return. +* If G=0, then CS=1 and SN=0. +* If F=0, then CS=0 and SN is chosen so that R is real. +* +* Arguments +* ========= +* +* F (input) COMPLEX*16 +* The first component of vector to be rotated. +* +* G (input) COMPLEX*16 +* The second component of vector to be rotated. +* +* CS (output) DOUBLE PRECISION +* The cosine of the rotation. +* +* SN (output) COMPLEX*16 +* The sine of the rotation. +* +* R (output) COMPLEX*16 +* The nonzero component of the rotated vector. +* +* Further Details +* ======= ======= +* +* 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel +* +* This version has a few statements commented out for thread safety +* (machine parameters are computed on each entry). 10 feb 03, SJH. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION TWO, ONE, ZERO + PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) + COMPLEX*16 CZERO + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. +* LOGICAL FIRST + INTEGER COUNT, I + DOUBLE PRECISION D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN, + $ SAFMN2, SAFMX2, SCALE + COMPLEX*16 FF, FS, GS +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLAPY2 + EXTERNAL DLAMCH, DLAPY2 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG, + $ MAX, SQRT +* .. +* .. Statement Functions .. + DOUBLE PRECISION ABS1, ABSSQ +* .. +* .. Save statement .. +* SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 +* .. +* .. Data statements .. +* DATA FIRST / .TRUE. / +* .. +* .. Statement Function definitions .. + ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) ) + ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2 +* .. +* .. Executable Statements .. +* +* IF( FIRST ) THEN + SAFMIN = DLAMCH( 'S' ) + EPS = DLAMCH( 'E' ) + SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / + $ LOG( DLAMCH( 'B' ) ) / TWO ) + SAFMX2 = ONE / SAFMN2 +* FIRST = .FALSE. +* END IF + SCALE = MAX( ABS1( F ), ABS1( G ) ) + FS = F + GS = G + COUNT = 0 + IF( SCALE.GE.SAFMX2 ) THEN + 10 CONTINUE + COUNT = COUNT + 1 + FS = FS*SAFMN2 + GS = GS*SAFMN2 + SCALE = SCALE*SAFMN2 + IF( SCALE.GE.SAFMX2 ) + $ GO TO 10 + ELSE IF( SCALE.LE.SAFMN2 ) THEN + IF( G.EQ.CZERO ) THEN + CS = ONE + SN = CZERO + R = F + RETURN + END IF + 20 CONTINUE + COUNT = COUNT - 1 + FS = FS*SAFMX2 + GS = GS*SAFMX2 + SCALE = SCALE*SAFMX2 + IF( SCALE.LE.SAFMN2 ) + $ GO TO 20 + END IF + F2 = ABSSQ( FS ) + G2 = ABSSQ( GS ) + IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN +* +* This is a rare case: F is very small. +* + IF( F.EQ.CZERO ) THEN + CS = ZERO + R = DLAPY2( DBLE( G ), DIMAG( G ) ) +* Do complex/real division explicitly with two real divisions + D = DLAPY2( DBLE( GS ), DIMAG( GS ) ) + SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D ) + RETURN + END IF + F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) ) +* G2 and G2S are accurate +* G2 is at least SAFMIN, and G2S is at least SAFMN2 + G2S = SQRT( G2 ) +* Error in CS from underflow in F2S is at most +* UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS +* If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN, +* and so CS .lt. sqrt(SAFMIN) +* If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN +* and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS) +* Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S + CS = F2S / G2S +* Make sure abs(FF) = 1 +* Do complex/real division explicitly with 2 real divisions + IF( ABS1( F ).GT.ONE ) THEN + D = DLAPY2( DBLE( F ), DIMAG( F ) ) + FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D ) + ELSE + DR = SAFMX2*DBLE( F ) + DI = SAFMX2*DIMAG( F ) + D = DLAPY2( DR, DI ) + FF = DCMPLX( DR / D, DI / D ) + END IF + SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S ) + R = CS*F + SN*G + ELSE +* +* This is the most common case. +* Neither F2 nor F2/G2 are less than SAFMIN +* F2S cannot overflow, and it is accurate +* + F2S = SQRT( ONE+G2 / F2 ) +* Do the F2S(real)*FS(complex) multiply with two real multiplies + R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) ) + CS = ONE / F2S + D = F2 + G2 +* Do complex/real division explicitly with two real divisions + SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D ) + SN = SN*DCONJG( GS ) + IF( COUNT.NE.0 ) THEN + IF( COUNT.GT.0 ) THEN + DO 30 I = 1, COUNT + R = R*SAFMX2 + 30 CONTINUE + ELSE + DO 40 I = 1, -COUNT + R = R*SAFMN2 + 40 CONTINUE + END IF + END IF + END IF + RETURN +* +* End of ZLARTG +* + END |