82
83
84
85
86
87
88 INTEGER INCX, N
89 COMPLEX*16 A
90
91
92 COMPLEX*16 X( * )
93
94
95
96
97
98 DOUBLE PRECISION ZERO, ONE
99 parameter( zero = 0.0d+0, one = 1.0d+0 )
100
101
102 DOUBLE PRECISION SAFMAX, SAFMIN, OV, AR, AI, ABSR, ABSI, UR, UI
103
104
105 DOUBLE PRECISION DLAMCH
106 COMPLEX*16 ZLADIV
108
109
111
112
113 INTRINSIC abs
114
115
116
117
118
119 IF( n.LE.0 )
120 $ RETURN
121
122
123
125 safmax = one / safmin
127
128
129
130 ar = dble( a )
131 ai = dimag( a )
132 absr = abs( ar )
133 absi = abs( ai )
134
135 IF( ai.EQ.zero ) THEN
136
137 CALL zdrscl( n, ar, x, incx )
138
139 ELSE IF( ar.EQ.zero ) THEN
140
141
142 IF( absi.GT.safmax ) THEN
143 CALL zdscal( n, safmin, x, incx )
144 CALL zscal( n, dcmplx( zero, -safmax / ai ), x, incx )
145 ELSE IF( absi.LT.safmin ) THEN
146 CALL zscal( n, dcmplx( zero, -safmin / ai ), x, incx )
147 CALL zdscal( n, safmax, x, incx )
148 ELSE
149 CALL zscal( n, dcmplx( zero, -one / ai ), x, incx )
150 END IF
151
152 ELSE
153
154
155
156
157
158
159
160 ur = ar + ai * ( ai / ar )
161 ui = ai + ar * ( ar / ai )
162
163 IF( (abs( ur ).LT.safmin).OR.(abs( ui ).LT.safmin) ) THEN
164
165 CALL zscal( n, dcmplx( safmin / ur, -safmin / ui ), x,
166 $ incx )
167 CALL zdscal( n, safmax, x, incx )
168 ELSE IF( (abs( ur ).GT.safmax).OR.(abs( ui ).GT.safmax) ) THEN
169 IF( (absr.GT.ov).OR.(absi.GT.ov) ) THEN
170
171 CALL zscal( n, dcmplx( one / ur, -one / ui ), x,
172 $ incx )
173 ELSE
174 CALL zdscal( n, safmin, x, incx )
175 IF( (abs( ur ).GT.ov).OR.(abs( ui ).GT.ov) ) THEN
176
177 IF( absr.GE.absi ) THEN
178
179 ur = (safmin * ar) + safmin * (ai * ( ai / ar ))
180 ui = (safmin * ai) + ar * ( (safmin * ar) / ai )
181 ELSE
182
183 ur = (safmin * ar) + ai * ( (safmin * ai) / ar )
184 ui = (safmin * ai) + safmin * (ar * ( ar / ai ))
185 END IF
186 CALL zscal( n, dcmplx( one / ur, -one / ui ), x,
187 $ incx )
188 ELSE
189 CALL zscal( n, dcmplx( safmax / ur, -safmax / ui ),
190 $ x, incx )
191 END IF
192 END IF
193 ELSE
194 CALL zscal( n, dcmplx( one / ur, -one / ui ), x, incx )
195 END IF
196 END IF
197
198 RETURN
199
200
201
complex *16 function zladiv(x, y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
double precision function dlamch(cmach)
DLAMCH
subroutine zdrscl(n, sa, sx, incx)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
subroutine dscal(n, da, dx, incx)
DSCAL
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine zscal(n, za, zx, incx)
ZSCAL