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