SUBROUTINE CLARTG( F, G, CS, SN, R ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * Jan 10, 2001 * * .. Scalar Arguments .. REAL CS COMPLEX F, G, R, SN * .. * * Purpose * ======= * * CLARTG 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 CROTG, except for * the following differences: * F and G are unchanged on return. * If F=0 and G=0, then CS=1, SN=0, and R=0. * If F .ne. 0 and G=0, then CS=1, SN=0, and R=F. * If F=0 and G .ne. 0, then CS=0, SN=conj(G)/abs(G), and R=abs(G). * If F .ne. 0 and G .ne. 0, then * CS = abs(F)/sqrt(F**2 + G**2) * SN = (F/abs(F))*conj(G)/sqrt(F**2 + G**2) * R = (F/abs(F))*sqrt(F**2 + G**2) * * If a NaN occurs in the input, then either real(R) and/or imag(R), * and possibly also CS, real(SN) and imag(SN), will be NaNs. * * If an infinity occurs in the input, then either real(R) and/or imag(R), * and possibly also CS, real(SN) and imag(SN), will be infinite or NaNs. * * The real routine SLARTG returns the same CS and SN (modulo roundoff) * if the inputs F and G are real. * * Arguments * ========= * * F (input) COMPLEX * The first component of vector to be rotated. * * G (input) COMPLEX * The second component of vector to be rotated. * * CS (output) REAL * The cosine of the rotation. * * SN (output) COMPLEX * The sine of the rotation. * * R (output) COMPLEX * The nonzero component of the rotated vector. * * ===================================================================== * * .. Parameters .. REAL FOUR, ONE, ZERO PARAMETER ( FOUR = 4.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) COMPLEX CZERO PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) * .. * .. Local Scalars .. LOGICAL FIRST, AGAIN INTEGER COUNT, I, MAXCNT, COUNTF, COUNTG REAL D1, EPS, F2, G2, SAFMIN, $ SAFMN2, SAFMX2, SAFMN4, SAFMX4, SAFMN, SAFMX, $ SCALEF, SCALEG, SCALEFG, FG2, SQREPS, ESFMN4, $ SCALE, SCALE2 COMPLEX FF, FS, GS * .. * .. External Functions .. REAL SLAMCH, SLAPY2 EXTERNAL SLAMCH, SLAPY2 * .. * .. Intrinsic Functions .. INTRINSIC ABS, AIMAG, CMPLX, CONJG, INT, LOG, MAX, REAL, $ SQRT * .. * .. Statement Functions .. REAL ABS1, ABSSQ * .. * .. Save statement .. SAVE FIRST, SAFMIN, EPS, SQREPS SAVE SAFMX2, SAFMX4, SAFMN2, SAFMN4, SAFMN, SAFMX * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Statement Function definitions .. ABS1( FF ) = MAX( ABS( REAL( FF ) ), ABS( AIMAG( FF ) ) ) ABSSQ( FF ) = REAL( FF )**2 + AIMAG( FF )**2 * .. * .. Executable Statements .. * IF( FIRST ) THEN * * On first call to CLARTG, compute * * SAFMN4 = (SAFMIN/EPS)**.25 rounded down to the nearest power * of the floating point radix * SAFMN2 = (SAFMIN/EPS)**.5 rounded down to the nearest power * of the floating point radix * * This means that scaling by SAFMN{2,4} and their * reciprocals SAFMX{2,4} causes no roundoff error * FIRST = .FALSE. SAFMIN = SLAMCH( 'S' ) EPS = SLAMCH( 'E' ) SQREPS = SQRT( EPS ) ESFMN4 = INT( LOG( SAFMIN / EPS ) / + LOG( SLAMCH( 'B' ) ) / FOUR ) SAFMN4 = SLAMCH( 'B' )**ESFMN4 SAFMN2 = SAFMN4**2 SAFMN = SAFMN2**2 SAFMX4 = ONE / SAFMN4 SAFMX2 = SAFMX4**2 SAFMX = SAFMX2**2 * * MAXCNT is the maximum number of times a nonzero number * can be scaled up/down by SAFMN4 before reaching 1 * MAXCNT = INT( -MAX( SLAMCH('L'), SLAMCH('N')+ONE-SLAMCH('M') ) + /ESFMN4 - ONE - ONE ) ENDIF * * If F (or G) contains a NaN, SCALEF (or SCALEG) might not SCALEF = ABS1( F ) SCALEG = ABS1( G ) IF( SCALEG.EQ.ZERO ) THEN * * Includes the case F=G=0 * Includes the case F is infinite or NaN: * If F is a NaN then R will be a NaN * If F is infinite then R will be infinite * or NaN if G is a NaN * May include cases where G is a NaN: * If G is a NaN then R will be a NaN * CS = ONE SN = CZERO * Add G to ensure that if G contains a NaN, so does R R = F + G ELSEIF( SCALEF.EQ.ZERO ) THEN * * G must be nonzero. * Includes the cases where G is infinite or NaN: * If G is a NaN then SN and R will be NaNs. * If G is infinite then SN will be a NAN and R will be infinite * May Include cases where F is a NaN: * If F is a NaN then R will be a NaN * CS = ZERO GS = G COUNT = 0 * No scaling if SCALEG is a NaN IF ( SCALEG .GT. SAFMX2 ) THEN 1 CONTINUE COUNT = COUNT + 1 GS = GS * SAFMN SCALEG = SCALEG * SAFMN * Keep scaling unless SCALEG is infinite IF ( SCALEG .GT. SAFMX2 .AND. COUNT .LE. MAXCNT ) GOTO 1 SCALE = SAFMX ELSEIF( SCALEG .LT. SAFMN2 ) THEN 2 CONTINUE COUNT = COUNT + 1 GS = GS * SAFMX SCALEG = SCALEG * SAFMX * Keep scaling unless SCALEG=0 because G contains a NaN IF ( SCALEG .LT. SAFMN2 .AND. COUNT .LE. MAXCNT ) GOTO 2 SCALE = SAFMN ENDIF D1 = SQRT( REAL(GS)**2 + AIMAG(GS)**2 ) R = D1 D1 = ONE/D1 SN = CMPLX( REAL(GS)*D1, -AIMAG(GS)*D1 ) DO 3 I = 1, COUNT R = CMPLX( REAL(R)*SCALE, ZERO ) 3 CONTINUE * Make sure that R contains a NaN if F does R = R + F ELSE * * Both F and G must be nonzero * IF( SCALEF.LE.SAFMX4 .AND. SCALEF.GE.SAFMN4 .AND. $ SCALEG.LE.SAFMX4 ) THEN * * Case 1: Neither F nor G too big or too small, minimal work * Neither F nor G can be infinite * If either F or G a NaN, then * CS, SR and R will be NaNs * F2 = ABSSQ(F) G2 = ABSSQ(G) FG2 = F2+G2 D1 = ONE/SQRT( F2*FG2 ) CS = F2*D1 FG2 = FG2 * D1 R = CMPLX( REAL(F)*FG2, AIMAG(F)*FG2 ) SN = CMPLX( REAL(F)*D1 , AIMAG(F)*D1 ) SN = CONJG(G) * SN ELSEIF( SCALEG .LT. SQREPS*SCALEF ) THEN * * Case 2: ABS(F)**2 + ABS(G)**2 rounds to ABS(F)**2 * F may be infinite but not G * Either F and/or G may be NaNs * CS = 1 always, even if F and/or G is a NaN * SN is a NaN if F or G is a NaN or infinite. * R is a NaN (or infinite) if F is a NaN (or infinite). * CS = ONE FS = F GS = G COUNTF = 0 COUNTG = 0 IF( SCALEF .GT. SAFMX2 ) THEN 10 CONTINUE COUNTF = COUNTF + 1 FS = FS * SAFMN SCALEF = SCALEF * SAFMN * Keep scaling unless SCALEF is infinite IF ( SCALEF .GT. SAFMX2 .AND. COUNTF .LE. MAXCNT) GOTO 10 ELSEIF( SCALEF .LT. SAFMN2 ) THEN 20 CONTINUE COUNTF = COUNTF - 1 FS = FS * SAFMX SCALEF = SCALEF * SAFMX * Keep scaling unless SCALEF=0 because F contains a NaN IF ( SCALEF .LT. SAFMN2 .AND. COUNTF.GE.-MAXCNT ) GOTO 20 ENDIF IF( SCALEG .GT. SAFMX2 ) THEN 30 CONTINUE COUNTG = COUNTG - 1 GS = GS * SAFMN SCALEG = SCALEG * SAFMN * SCALEG finite so scaling must terminate IF ( SCALEG .GT. SAFMX2 ) GOTO 30 ELSEIF( SCALEG .LT. SAFMN2 ) THEN 40 CONTINUE COUNTG = COUNTG + 1 GS = GS * SAFMX SCALEG = SCALEG * SAFMX * Keep scaling unless SCALEG=0 because G contains a NaN IF ( SCALEG .LT. SAFMN2 .AND. COUNTG .LE.MAXCNT ) GOTO 40 ENDIF D1 = ONE/(REAL(FS)**2 + AIMAG(FS)**2) SN = FS * CONJG(GS) * SN will be a NaN if F is infinite or a NaN SN = CMPLX( REAL(SN)*D1 , AIMAG(SN)*D1 ) COUNT = COUNTF + COUNTG IF( COUNT .GT. 0 ) THEN DO 50 I = 1, COUNT SN = CMPLX( REAL(SN)*SAFMN , AIMAG(SN)*SAFMN ) 50 CONTINUE ELSEIF( COUNT .LT. 0 ) THEN DO 60 I = 1, -COUNT SN = CMPLX( REAL(SN)*SAFMX , AIMAG(SN)*SAFMX ) 60 CONTINUE ENDIF * Make sure R contains a NaN if G does R = F + SN*G ELSEIF( SCALEF .LT. SQREPS*SCALEG ) THEN * * Case 3: ABS(F)**2 + ABS(G)**2 rounds to ABS(G)**2 * G may be infinite but not F * Either F or G may be a NaN, in which case * CS, SN and R are NaNs * SN and R are NaNs if G is infinite * FS = F GS = G COUNTF = 0 COUNTG = 0 IF( SCALEF .GT. SAFMX4 ) THEN 70 CONTINUE COUNTF = COUNTF + 1 FS = FS * SAFMN2 SCALEF = SCALEF * SAFMN2 * SCALEF finite so scaling must terminate IF ( SCALEF .GT. SAFMX4 ) GOTO 70 ELSEIF( SCALEF .LT. SAFMN4 ) THEN 80 CONTINUE COUNTF = COUNTF - 1 FS = FS * SAFMX2 SCALEF = SCALEF * SAFMX2 * Keep scaling unless SCALEF=0 because F contains a NaN IF ( SCALEF .LT. SAFMN4 .AND. COUNTF .GE.-MAXCNT) GOTO 80 ENDIF IF( SCALEG .GT. SAFMX4 ) THEN 90 CONTINUE COUNTG = COUNTG + 1 GS = GS * SAFMN2 SCALEG = SCALEG * SAFMN2 * Keep scaling unless SCALEG is infinite IF ( SCALEG .GT. SAFMX4 .AND. COUNTG .LE. MAXCNT) GOTO 90 ELSEIF( SCALEG .LT. SAFMN4 ) THEN 100 CONTINUE COUNTG = COUNTG - 1 GS = GS * SAFMX2 SCALEG = SCALEG * SAFMX2 * SCALEG cannot be zero so scaling must terminate IF ( SCALEG .LT. SAFMN4 ) GOTO 100 ENDIF F2 = REAL(FS)**2 + AIMAG(FS)**2 G2 = REAL(GS)**2 + AIMAG(GS)**2 D1 = ONE/SQRT( F2*G2 ) CS = F2*D1 SN = FS * CONJG(GS) * SN will be a NaN if G is infinite SN = CMPLX( REAL(SN)*D1 , AIMAG(SN)*D1 ) D1 = G2*D1 * R will be a NaN if G is infinite R = CMPLX( REAL(FS)*D1 , AIMAG(FS)*D1 ) COUNT = COUNTF - COUNTG IF( COUNT .GT. 0 ) THEN DO 110 I = 1, COUNT CS = CS*SAFMX2 110 CONTINUE ELSEIF( COUNT .LT. 0 ) THEN DO 120 I = 1, -COUNT CS = CS*SAFMN2 120 CONTINUE ENDIF IF( COUNTG .GT. 0 ) THEN DO 130 I = 1, COUNTG R = CMPLX( REAL(R)*SAFMX2, AIMAG(R)*SAFMX2 ) 130 CONTINUE ELSEIF( COUNTG .LT. 0 ) THEN DO 140 I = 1, -COUNTG R = CMPLX( REAL(R)*SAFMN2, AIMAG(R)*SAFMN2 ) 140 CONTINUE ENDIF ELSE * * Case 4: Scale F and G up or down and use formula from Case 1 * Both F and G may be simultaneously infinite, * in which case CS, SN and R are all NaNs * Either F or G may be a NaN, in which case * CS, SN and R are all NaNs * FS = F GS = G COUNT = 0 AGAIN = .FALSE. SCALEFG = MAX( SCALEF, SCALEG ) IF( SCALEFG .GT. ONE) THEN * IF (SCALEFG .LE. SAFMX2) THEN SCALE2 = SAFMX4 AGAIN = .TRUE. FS = FS * SAFMN4 GS = GS * SAFMN4 GOTO 153 ENDIF SCALE = SAFMX2 FS = FS * SAFMN2 GS = GS * SAFMN2 SCALEFG = SCALEFG * SAFMN2 COUNT = COUNT + 1 * * Keep scaling unless SCALEFG is infinite 151 CONTINUE IF( SCALEFG .LE. SAFMX2 .OR. COUNT .GT. MAXCNT) GOTO 152 FS = FS * SAFMN2 GS = GS * SAFMN2 SCALEFG = SCALEFG * SAFMN2 COUNT = COUNT + 1 GOTO 151 152 CONTINUE IF( SCALEFG .GT. SAFMX4) THEN SCALE2 = SAFMX4 AGAIN = .TRUE. FS = FS * SAFMN4 GS = GS * SAFMN4 ENDIF ELSE * SCALEFG might be a NaN * IF (SCALEF .GE. SAFMN2) THEN SCALE2 = SAFMN4 AGAIN = .TRUE. FS = FS * SAFMX4 GS = GS * SAFMX4 GOTO 153 ENDIF SCALE = SAFMN2 FS = FS * SAFMX2 GS = GS * SAFMX2 SCALEF = SCALEF * SAFMX2 COUNT = COUNT + 1 * 160 CONTINUE IF( SCALEF .GE. SAFMN2 .OR. COUNT .GT. MAXCNT ) GOTO 162 FS = FS * SAFMX2 GS = GS * SAFMX2 SCALEF = SCALEF * SAFMX2 COUNT = COUNT + 1 162 CONTINUE IF( SCALEF .LT. SAFMN4 ) THEN SCALE2 = SAFMN4 AGAIN = .TRUE. FS = FS * SAFMX4 GS = GS * SAFMX4 ENDIF ENDIF * 153 CONTINUE F2 = ABSSQ(FS) G2 = ABSSQ(GS) FG2 = F2+G2 D1 = ONE/SQRT( F2*FG2 ) * CS will be a NaN if both F and G infinite CS = F2*D1 FG2 = FG2 * D1 R = CMPLX( REAL(FS)*FG2, AIMAG(FS)*FG2 ) SN = CMPLX( REAL(FS)*D1 , AIMAG(FS)*D1 ) * SN will be a NaN if both F and G infinite SN = CONJG(GS) * SN DO 170 I = 1, COUNT R = CMPLX( REAL(R) * SCALE, AIMAG(R) * SCALE ) 170 CONTINUE IF ( AGAIN ) $ R = CMPLX( REAL(R) * SCALE2, AIMAG(R) * SCALE2 ) ENDIF ENDIF RETURN * * End of CLARTG * END