104
105
106
107
108
109
110 INTEGER INCX, N
111 COMPLEX ALPHA, TAU
112
113
114 COMPLEX X( * )
115
116
117
118
119
120 REAL TWO, ONE, ZERO
121 parameter( two = 2.0e+0, one = 1.0e+0, zero = 0.0e+0 )
122
123
124 INTEGER J, KNT
125 REAL ALPHI, ALPHR, BETA, BIGNUM, EPS, SMLNUM, XNORM
126 COMPLEX SAVEALPHA
127
128
129 REAL SCNRM2, SLAMCH, SLAPY3, SLAPY2
130 COMPLEX CLADIV
132
133
134 INTRINSIC abs, aimag, cmplx, real, sign
135
136
138
139
140
141 IF( n.LE.0 ) THEN
142 tau = zero
143 RETURN
144 END IF
145
146 eps =
slamch(
'Precision' )
147 xnorm =
scnrm2( n-1, x, incx )
148 alphr = real( alpha )
149 alphi = aimag( alpha )
150
151 IF( xnorm.LE.eps*abs(alpha) ) THEN
152
153
154
155 IF( alphi.EQ.zero ) THEN
156 IF( alphr.GE.zero ) THEN
157
158
159
160 tau = zero
161 ELSE
162
163
164 tau = two
165 DO j = 1, n-1
166 x( 1 + (j-1)*incx ) = zero
167 END DO
168 alpha = -alpha
169 END IF
170 ELSE
171
172 xnorm =
slapy2( alphr, alphi )
173 tau = cmplx( one - alphr / xnorm, -alphi / xnorm )
174 DO j = 1, n-1
175 x( 1 + (j-1)*incx ) = zero
176 END DO
177 alpha = xnorm
178 END IF
179 ELSE
180
181
182
183 beta = sign(
slapy3( alphr, alphi, xnorm ), alphr )
185 bignum = one / smlnum
186
187 knt = 0
188 IF( abs( beta ).LT.smlnum ) THEN
189
190
191
192 10 CONTINUE
193 knt = knt + 1
194 CALL csscal( n-1, bignum, x, incx )
195 beta = beta*bignum
196 alphi = alphi*bignum
197 alphr = alphr*bignum
198 IF( (abs( beta ).LT.smlnum) .AND. (knt .LT. 20) )
199 $ GO TO 10
200
201
202
203 xnorm =
scnrm2( n-1, x, incx )
204 alpha = cmplx( alphr, alphi )
205 beta = sign(
slapy3( alphr, alphi, xnorm ), alphr )
206 END IF
207 savealpha = alpha
208 alpha = alpha + beta
209 IF( beta.LT.zero ) THEN
210 beta = -beta
211 tau = -alpha / beta
212 ELSE
213 alphr = alphi * (alphi/real( alpha ))
214 alphr = alphr + xnorm * (xnorm/real( alpha ))
215 tau = cmplx( alphr/beta, -alphi/beta )
216 alpha = cmplx( -alphr, alphi )
217 END IF
218 alpha =
cladiv( cmplx( one ), alpha )
219
220 IF ( abs(tau).LE.smlnum ) THEN
221
222
223
224
225
226
227
228
229 alphr = real( savealpha )
230 alphi = aimag( savealpha )
231 IF( alphi.EQ.zero ) THEN
232 IF( alphr.GE.zero ) THEN
233 tau = zero
234 ELSE
235 tau = two
236 DO j = 1, n-1
237 x( 1 + (j-1)*incx ) = zero
238 END DO
239 beta = real( -savealpha )
240 END IF
241 ELSE
242 xnorm =
slapy2( alphr, alphi )
243 tau = cmplx( one - alphr / xnorm, -alphi / xnorm )
244 DO j = 1, n-1
245 x( 1 + (j-1)*incx ) = zero
246 END DO
247 beta = xnorm
248 END IF
249
250 ELSE
251
252
253
254 CALL cscal( n-1, alpha, x, incx )
255
256 END IF
257
258
259
260 DO 20 j = 1, knt
261 beta = beta*smlnum
262 20 CONTINUE
263 alpha = beta
264 END IF
265
266 RETURN
267
268
269
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
real function slamch(cmach)
SLAMCH
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
real function slapy3(x, y, z)
SLAPY3 returns sqrt(x2+y2+z2).
real(wp) function scnrm2(n, x, incx)
SCNRM2
subroutine csscal(n, sa, cx, incx)
CSSCAL
subroutine cscal(n, ca, cx, incx)
CSCAL