114 COMPLEX*16 alpha, tau
123 DOUBLE PRECISION two, one, zero
124 parameter( two = 2.0d+0, one = 1.0d+0, zero = 0.0d+0 )
128 DOUBLE PRECISION alphi, alphr, beta, bignum, smlnum, xnorm
137 INTRINSIC abs, dble, dcmplx, dimag, sign
150 alphr = dble( alpha )
151 alphi = dimag( alpha )
153 IF( xnorm.EQ.zero )
THEN
157 IF( alphi.EQ.zero )
THEN
158 IF( alphr.GE.zero )
THEN
168 x( 1 + (
j-1)*incx ) = zero
174 xnorm =
dlapy2( alphr, alphi )
175 tau = dcmplx( one - alphr / xnorm, -alphi / xnorm )
177 x( 1 + (
j-1)*incx ) = zero
185 beta = sign(
dlapy3( alphr, alphi, xnorm ), alphr )
187 bignum = one / smlnum
190 IF( abs( beta ).LT.smlnum )
THEN
196 CALL
zdscal( n-1, bignum,
x, incx )
200 IF( abs( beta ).LT.smlnum )
206 alpha = dcmplx( alphr, alphi )
207 beta = sign(
dlapy3( alphr, alphi, xnorm ), alphr )
211 IF( beta.LT.zero )
THEN
215 alphr = alphi * (alphi/dble( alpha ))
216 alphr = alphr + xnorm * (xnorm/dble( alpha ))
217 tau = dcmplx( alphr/beta, -alphi/beta )
218 alpha = dcmplx( -alphr, alphi )
220 alpha =
zladiv( dcmplx( one ), alpha )
222 IF ( abs(tau).LE.smlnum )
THEN
231 alphr = dble( savealpha )
232 alphi = dimag( savealpha )
233 IF( alphi.EQ.zero )
THEN
234 IF( alphr.GE.zero )
THEN
239 x( 1 + (
j-1)*incx ) = zero
244 xnorm =
dlapy2( alphr, alphi )
245 tau = dcmplx( one - alphr / xnorm, -alphi / xnorm )
247 x( 1 + (
j-1)*incx ) = zero
256 CALL
zscal( n-1, alpha,
x, incx )
COMPLEX *16 function zladiv(X, Y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
DOUBLE PRECISION function dlapy3(X, Y, Z)
DLAPY3 returns sqrt(x2+y2+z2).
subroutine zlarfgp(N, ALPHA, X, INCX, TAU)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negatibe beta.
DOUBLE PRECISION function dznrm2(N, X, INCX)
DZNRM2
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
DOUBLE PRECISION function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL