```001:       SUBROUTINE ZLARTG( F, G, CS, SN, R )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       DOUBLE PRECISION   CS
009:       COMPLEX*16         F, G, R, SN
010: *     ..
011: *
012: *  Purpose
013: *  =======
014: *
015: *  ZLARTG generates a plane rotation so that
016: *
017: *     [  CS  SN  ]     [ F ]     [ R ]
018: *     [  __      ]  .  [   ]  =  [   ]   where CS**2 + |SN|**2 = 1.
019: *     [ -SN  CS  ]     [ G ]     [ 0 ]
020: *
021: *  This is a faster version of the BLAS1 routine ZROTG, except for
022: *  the following differences:
023: *     F and G are unchanged on return.
024: *     If G=0, then CS=1 and SN=0.
025: *     If F=0, then CS=0 and SN is chosen so that R is real.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  F       (input) COMPLEX*16
031: *          The first component of vector to be rotated.
032: *
033: *  G       (input) COMPLEX*16
034: *          The second component of vector to be rotated.
035: *
036: *  CS      (output) DOUBLE PRECISION
037: *          The cosine of the rotation.
038: *
039: *  SN      (output) COMPLEX*16
040: *          The sine of the rotation.
041: *
042: *  R       (output) COMPLEX*16
043: *          The nonzero component of the rotated vector.
044: *
045: *  Further Details
046: *  ======= =======
047: *
048: *  3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel
049: *
050: *  This version has a few statements commented out for thread safety
051: *  (machine parameters are computed on each entry). 10 feb 03, SJH.
052: *
053: *  =====================================================================
054: *
055: *     .. Parameters ..
056:       DOUBLE PRECISION   TWO, ONE, ZERO
057:       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
058:       COMPLEX*16         CZERO
059:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
060: *     ..
061: *     .. Local Scalars ..
062: *     LOGICAL            FIRST
063:       INTEGER            COUNT, I
064:       DOUBLE PRECISION   D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
065:      \$                   SAFMN2, SAFMX2, SCALE
066:       COMPLEX*16         FF, FS, GS
067: *     ..
068: *     .. External Functions ..
069:       DOUBLE PRECISION   DLAMCH, DLAPY2
070:       EXTERNAL           DLAMCH, DLAPY2
071: *     ..
072: *     .. Intrinsic Functions ..
073:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
074:      \$                   MAX, SQRT
075: *     ..
076: *     .. Statement Functions ..
077:       DOUBLE PRECISION   ABS1, ABSSQ
078: *     ..
079: *     .. Save statement ..
080: *     SAVE               FIRST, SAFMX2, SAFMIN, SAFMN2
081: *     ..
082: *     .. Data statements ..
083: *     DATA               FIRST / .TRUE. /
084: *     ..
085: *     .. Statement Function definitions ..
086:       ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
087:       ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
088: *     ..
089: *     .. Executable Statements ..
090: *
091: *     IF( FIRST ) THEN
092:          SAFMIN = DLAMCH( 'S' )
093:          EPS = DLAMCH( 'E' )
094:          SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
095:      \$            LOG( DLAMCH( 'B' ) ) / TWO )
096:          SAFMX2 = ONE / SAFMN2
097: *        FIRST = .FALSE.
098: *     END IF
099:       SCALE = MAX( ABS1( F ), ABS1( G ) )
100:       FS = F
101:       GS = G
102:       COUNT = 0
103:       IF( SCALE.GE.SAFMX2 ) THEN
104:    10    CONTINUE
105:          COUNT = COUNT + 1
106:          FS = FS*SAFMN2
107:          GS = GS*SAFMN2
108:          SCALE = SCALE*SAFMN2
109:          IF( SCALE.GE.SAFMX2 )
110:      \$      GO TO 10
111:       ELSE IF( SCALE.LE.SAFMN2 ) THEN
112:          IF( G.EQ.CZERO ) THEN
113:             CS = ONE
114:             SN = CZERO
115:             R = F
116:             RETURN
117:          END IF
118:    20    CONTINUE
119:          COUNT = COUNT - 1
120:          FS = FS*SAFMX2
121:          GS = GS*SAFMX2
122:          SCALE = SCALE*SAFMX2
123:          IF( SCALE.LE.SAFMN2 )
124:      \$      GO TO 20
125:       END IF
126:       F2 = ABSSQ( FS )
127:       G2 = ABSSQ( GS )
128:       IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
129: *
130: *        This is a rare case: F is very small.
131: *
132:          IF( F.EQ.CZERO ) THEN
133:             CS = ZERO
134:             R = DLAPY2( DBLE( G ), DIMAG( G ) )
135: *           Do complex/real division explicitly with two real divisions
136:             D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
137:             SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
138:             RETURN
139:          END IF
140:          F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
141: *        G2 and G2S are accurate
142: *        G2 is at least SAFMIN, and G2S is at least SAFMN2
143:          G2S = SQRT( G2 )
144: *        Error in CS from underflow in F2S is at most
145: *        UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
146: *        If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
147: *        and so CS .lt. sqrt(SAFMIN)
148: *        If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
149: *        and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
150: *        Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
151:          CS = F2S / G2S
152: *        Make sure abs(FF) = 1
153: *        Do complex/real division explicitly with 2 real divisions
154:          IF( ABS1( F ).GT.ONE ) THEN
155:             D = DLAPY2( DBLE( F ), DIMAG( F ) )
156:             FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
157:          ELSE
158:             DR = SAFMX2*DBLE( F )
159:             DI = SAFMX2*DIMAG( F )
160:             D = DLAPY2( DR, DI )
161:             FF = DCMPLX( DR / D, DI / D )
162:          END IF
163:          SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
164:          R = CS*F + SN*G
165:       ELSE
166: *
167: *        This is the most common case.
168: *        Neither F2 nor F2/G2 are less than SAFMIN
169: *        F2S cannot overflow, and it is accurate
170: *
171:          F2S = SQRT( ONE+G2 / F2 )
172: *        Do the F2S(real)*FS(complex) multiply with two real multiplies
173:          R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
174:          CS = ONE / F2S
175:          D = F2 + G2
176: *        Do complex/real division explicitly with two real divisions
177:          SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
178:          SN = SN*DCONJG( GS )
179:          IF( COUNT.NE.0 ) THEN
180:             IF( COUNT.GT.0 ) THEN
181:                DO 30 I = 1, COUNT
182:                   R = R*SAFMX2
183:    30          CONTINUE
184:             ELSE
185:                DO 40 I = 1, -COUNT
186:                   R = R*SAFMN2
187:    40          CONTINUE
188:             END IF
189:          END IF
190:       END IF
191:       RETURN
192: *
193: *     End of ZLARTG
194: *
195:       END
196: ```