1 SUBROUTINE ZLARGV( N, X, INCX, Y, INCY, C, INCC )
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INCC, INCX, INCY, N
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION C( * )
13 COMPLEX*16 X( * ), Y( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLARGV generates a vector of complex plane rotations with real
20 * cosines, determined by elements of the complex vectors x and y.
21 * For i = 1,2,...,n
22 *
23 * ( c(i) s(i) ) ( x(i) ) = ( r(i) )
24 * ( -conjg(s(i)) c(i) ) ( y(i) ) = ( 0 )
25 *
26 * where c(i)**2 + ABS(s(i))**2 = 1
27 *
28 * The following conventions are used (these are the same as in ZLARTG,
29 * but differ from the BLAS1 routine ZROTG):
30 * If y(i)=0, then c(i)=1 and s(i)=0.
31 * If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.
32 *
33 * Arguments
34 * =========
35 *
36 * N (input) INTEGER
37 * The number of plane rotations to be generated.
38 *
39 * X (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCX)
40 * On entry, the vector x.
41 * On exit, x(i) is overwritten by r(i), for i = 1,...,n.
42 *
43 * INCX (input) INTEGER
44 * The increment between elements of X. INCX > 0.
45 *
46 * Y (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCY)
47 * On entry, the vector y.
48 * On exit, the sines of the plane rotations.
49 *
50 * INCY (input) INTEGER
51 * The increment between elements of Y. INCY > 0.
52 *
53 * C (output) DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
54 * The cosines of the plane rotations.
55 *
56 * INCC (input) INTEGER
57 * The increment between elements of C. INCC > 0.
58 *
59 * Further Details
60 * ======= =======
61 *
62 * 6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel
63 *
64 * This version has a few statements commented out for thread safety
65 * (machine parameters are computed on each entry). 10 feb 03, SJH.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 DOUBLE PRECISION TWO, ONE, ZERO
71 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
72 COMPLEX*16 CZERO
73 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
74 * ..
75 * .. Local Scalars ..
76 * LOGICAL FIRST
77
78 INTEGER COUNT, I, IC, IX, IY, J
79 DOUBLE PRECISION CS, D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
80 $ SAFMN2, SAFMX2, SCALE
81 COMPLEX*16 F, FF, FS, G, GS, R, SN
82 * ..
83 * .. External Functions ..
84 DOUBLE PRECISION DLAMCH, DLAPY2
85 EXTERNAL DLAMCH, DLAPY2
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
89 $ MAX, SQRT
90 * ..
91 * .. Statement Functions ..
92 DOUBLE PRECISION ABS1, ABSSQ
93 * ..
94 * .. Save statement ..
95 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
96 * ..
97 * .. Data statements ..
98 * DATA FIRST / .TRUE. /
99 * ..
100 * .. Statement Function definitions ..
101 ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
102 ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
103 * ..
104 * .. Executable Statements ..
105 *
106 * IF( FIRST ) THEN
107 * FIRST = .FALSE.
108 SAFMIN = DLAMCH( 'S' )
109 EPS = DLAMCH( 'E' )
110 SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
111 $ LOG( DLAMCH( 'B' ) ) / TWO )
112 SAFMX2 = ONE / SAFMN2
113 * END IF
114 IX = 1
115 IY = 1
116 IC = 1
117 DO 60 I = 1, N
118 F = X( IX )
119 G = Y( IY )
120 *
121 * Use identical algorithm as in ZLARTG
122 *
123 SCALE = MAX( ABS1( F ), ABS1( G ) )
124 FS = F
125 GS = G
126 COUNT = 0
127 IF( SCALE.GE.SAFMX2 ) THEN
128 10 CONTINUE
129 COUNT = COUNT + 1
130 FS = FS*SAFMN2
131 GS = GS*SAFMN2
132 SCALE = SCALE*SAFMN2
133 IF( SCALE.GE.SAFMX2 )
134 $ GO TO 10
135 ELSE IF( SCALE.LE.SAFMN2 ) THEN
136 IF( G.EQ.CZERO ) THEN
137 CS = ONE
138 SN = CZERO
139 R = F
140 GO TO 50
141 END IF
142 20 CONTINUE
143 COUNT = COUNT - 1
144 FS = FS*SAFMX2
145 GS = GS*SAFMX2
146 SCALE = SCALE*SAFMX2
147 IF( SCALE.LE.SAFMN2 )
148 $ GO TO 20
149 END IF
150 F2 = ABSSQ( FS )
151 G2 = ABSSQ( GS )
152 IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
153 *
154 * This is a rare case: F is very small.
155 *
156 IF( F.EQ.CZERO ) THEN
157 CS = ZERO
158 R = DLAPY2( DBLE( G ), DIMAG( G ) )
159 * Do complex/real division explicitly with two real
160 * divisions
161 D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
162 SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
163 GO TO 50
164 END IF
165 F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
166 * G2 and G2S are accurate
167 * G2 is at least SAFMIN, and G2S is at least SAFMN2
168 G2S = SQRT( G2 )
169 * Error in CS from underflow in F2S is at most
170 * UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
171 * If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
172 * and so CS .lt. sqrt(SAFMIN)
173 * If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
174 * and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
175 * Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
176 CS = F2S / G2S
177 * Make sure abs(FF) = 1
178 * Do complex/real division explicitly with 2 real divisions
179 IF( ABS1( F ).GT.ONE ) THEN
180 D = DLAPY2( DBLE( F ), DIMAG( F ) )
181 FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
182 ELSE
183 DR = SAFMX2*DBLE( F )
184 DI = SAFMX2*DIMAG( F )
185 D = DLAPY2( DR, DI )
186 FF = DCMPLX( DR / D, DI / D )
187 END IF
188 SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
189 R = CS*F + SN*G
190 ELSE
191 *
192 * This is the most common case.
193 * Neither F2 nor F2/G2 are less than SAFMIN
194 * F2S cannot overflow, and it is accurate
195 *
196 F2S = SQRT( ONE+G2 / F2 )
197 * Do the F2S(real)*FS(complex) multiply with two real
198 * multiplies
199 R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
200 CS = ONE / F2S
201 D = F2 + G2
202 * Do complex/real division explicitly with two real divisions
203 SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
204 SN = SN*DCONJG( GS )
205 IF( COUNT.NE.0 ) THEN
206 IF( COUNT.GT.0 ) THEN
207 DO 30 J = 1, COUNT
208 R = R*SAFMX2
209 30 CONTINUE
210 ELSE
211 DO 40 J = 1, -COUNT
212 R = R*SAFMN2
213 40 CONTINUE
214 END IF
215 END IF
216 END IF
217 50 CONTINUE
218 C( IC ) = CS
219 Y( IY ) = SN
220 X( IX ) = R
221 IC = IC + INCC
222 IY = IY + INCY
223 IX = IX + INCX
224 60 CONTINUE
225 RETURN
226 *
227 * End of ZLARGV
228 *
229 END
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INCC, INCX, INCY, N
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION C( * )
13 COMPLEX*16 X( * ), Y( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLARGV generates a vector of complex plane rotations with real
20 * cosines, determined by elements of the complex vectors x and y.
21 * For i = 1,2,...,n
22 *
23 * ( c(i) s(i) ) ( x(i) ) = ( r(i) )
24 * ( -conjg(s(i)) c(i) ) ( y(i) ) = ( 0 )
25 *
26 * where c(i)**2 + ABS(s(i))**2 = 1
27 *
28 * The following conventions are used (these are the same as in ZLARTG,
29 * but differ from the BLAS1 routine ZROTG):
30 * If y(i)=0, then c(i)=1 and s(i)=0.
31 * If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.
32 *
33 * Arguments
34 * =========
35 *
36 * N (input) INTEGER
37 * The number of plane rotations to be generated.
38 *
39 * X (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCX)
40 * On entry, the vector x.
41 * On exit, x(i) is overwritten by r(i), for i = 1,...,n.
42 *
43 * INCX (input) INTEGER
44 * The increment between elements of X. INCX > 0.
45 *
46 * Y (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCY)
47 * On entry, the vector y.
48 * On exit, the sines of the plane rotations.
49 *
50 * INCY (input) INTEGER
51 * The increment between elements of Y. INCY > 0.
52 *
53 * C (output) DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
54 * The cosines of the plane rotations.
55 *
56 * INCC (input) INTEGER
57 * The increment between elements of C. INCC > 0.
58 *
59 * Further Details
60 * ======= =======
61 *
62 * 6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel
63 *
64 * This version has a few statements commented out for thread safety
65 * (machine parameters are computed on each entry). 10 feb 03, SJH.
66 *
67 * =====================================================================
68 *
69 * .. Parameters ..
70 DOUBLE PRECISION TWO, ONE, ZERO
71 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
72 COMPLEX*16 CZERO
73 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
74 * ..
75 * .. Local Scalars ..
76 * LOGICAL FIRST
77
78 INTEGER COUNT, I, IC, IX, IY, J
79 DOUBLE PRECISION CS, D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
80 $ SAFMN2, SAFMX2, SCALE
81 COMPLEX*16 F, FF, FS, G, GS, R, SN
82 * ..
83 * .. External Functions ..
84 DOUBLE PRECISION DLAMCH, DLAPY2
85 EXTERNAL DLAMCH, DLAPY2
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
89 $ MAX, SQRT
90 * ..
91 * .. Statement Functions ..
92 DOUBLE PRECISION ABS1, ABSSQ
93 * ..
94 * .. Save statement ..
95 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
96 * ..
97 * .. Data statements ..
98 * DATA FIRST / .TRUE. /
99 * ..
100 * .. Statement Function definitions ..
101 ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
102 ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
103 * ..
104 * .. Executable Statements ..
105 *
106 * IF( FIRST ) THEN
107 * FIRST = .FALSE.
108 SAFMIN = DLAMCH( 'S' )
109 EPS = DLAMCH( 'E' )
110 SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
111 $ LOG( DLAMCH( 'B' ) ) / TWO )
112 SAFMX2 = ONE / SAFMN2
113 * END IF
114 IX = 1
115 IY = 1
116 IC = 1
117 DO 60 I = 1, N
118 F = X( IX )
119 G = Y( IY )
120 *
121 * Use identical algorithm as in ZLARTG
122 *
123 SCALE = MAX( ABS1( F ), ABS1( G ) )
124 FS = F
125 GS = G
126 COUNT = 0
127 IF( SCALE.GE.SAFMX2 ) THEN
128 10 CONTINUE
129 COUNT = COUNT + 1
130 FS = FS*SAFMN2
131 GS = GS*SAFMN2
132 SCALE = SCALE*SAFMN2
133 IF( SCALE.GE.SAFMX2 )
134 $ GO TO 10
135 ELSE IF( SCALE.LE.SAFMN2 ) THEN
136 IF( G.EQ.CZERO ) THEN
137 CS = ONE
138 SN = CZERO
139 R = F
140 GO TO 50
141 END IF
142 20 CONTINUE
143 COUNT = COUNT - 1
144 FS = FS*SAFMX2
145 GS = GS*SAFMX2
146 SCALE = SCALE*SAFMX2
147 IF( SCALE.LE.SAFMN2 )
148 $ GO TO 20
149 END IF
150 F2 = ABSSQ( FS )
151 G2 = ABSSQ( GS )
152 IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
153 *
154 * This is a rare case: F is very small.
155 *
156 IF( F.EQ.CZERO ) THEN
157 CS = ZERO
158 R = DLAPY2( DBLE( G ), DIMAG( G ) )
159 * Do complex/real division explicitly with two real
160 * divisions
161 D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
162 SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
163 GO TO 50
164 END IF
165 F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
166 * G2 and G2S are accurate
167 * G2 is at least SAFMIN, and G2S is at least SAFMN2
168 G2S = SQRT( G2 )
169 * Error in CS from underflow in F2S is at most
170 * UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
171 * If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
172 * and so CS .lt. sqrt(SAFMIN)
173 * If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
174 * and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
175 * Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
176 CS = F2S / G2S
177 * Make sure abs(FF) = 1
178 * Do complex/real division explicitly with 2 real divisions
179 IF( ABS1( F ).GT.ONE ) THEN
180 D = DLAPY2( DBLE( F ), DIMAG( F ) )
181 FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
182 ELSE
183 DR = SAFMX2*DBLE( F )
184 DI = SAFMX2*DIMAG( F )
185 D = DLAPY2( DR, DI )
186 FF = DCMPLX( DR / D, DI / D )
187 END IF
188 SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
189 R = CS*F + SN*G
190 ELSE
191 *
192 * This is the most common case.
193 * Neither F2 nor F2/G2 are less than SAFMIN
194 * F2S cannot overflow, and it is accurate
195 *
196 F2S = SQRT( ONE+G2 / F2 )
197 * Do the F2S(real)*FS(complex) multiply with two real
198 * multiplies
199 R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
200 CS = ONE / F2S
201 D = F2 + G2
202 * Do complex/real division explicitly with two real divisions
203 SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
204 SN = SN*DCONJG( GS )
205 IF( COUNT.NE.0 ) THEN
206 IF( COUNT.GT.0 ) THEN
207 DO 30 J = 1, COUNT
208 R = R*SAFMX2
209 30 CONTINUE
210 ELSE
211 DO 40 J = 1, -COUNT
212 R = R*SAFMN2
213 40 CONTINUE
214 END IF
215 END IF
216 END IF
217 50 CONTINUE
218 C( IC ) = CS
219 Y( IY ) = SN
220 X( IX ) = R
221 IC = IC + INCC
222 IY = IY + INCY
223 IX = IX + INCX
224 60 CONTINUE
225 RETURN
226 *
227 * End of ZLARGV
228 *
229 END